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1 SUMMARY 

During the period of this support, a new control system design and analysis method has been 
studied. This approach deals with control systems containing uncertainties that are represented in 
terms of its transfer function parameters. Such a representation of the control system is common and 
many physical parameter variations fall into this type of uncertainty. Techniques developed here are 
capable of providing nonconservative analysis of such control systems with parameter variations. We 
have also developed techniques to deal with control systems when their state space representations 
are given rather than transfer functions. In this case, the plant parameters will appear as entries of 
state space matrices. Finally, a system modeling technique to construct such systems from the raw 
input - output frequency domain data has been developed. 
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3 INTRODUCTION 


The subject of robust control began to receive worldwide attention since the late 1970 ’s. In particu- 
lar, since the seminal theorem of Kharitonov, the problem of robustness under parametric uncertainty 
became one of the major subjects in which researchers heavily invested their efforts. Through this 
project, we aimed to achieve meaningful analysis and design mechanisms for a system whose un- 
certainties are modeled as parameters of the transfer function coefficients. At present the theory 
is capable of determining the largest parametric stability margin in various forms of measures, the 
worst case gain and phase margins, the worst case unstructured (ifoo) stability margin, and the 
nonlinear sector bounded stability margin while the given parametric uncertainty is simultaneously 
presented. These novel analysis techniques are equally applicable to systems whose transfer function 
coefficients are linear or multilinear functions of parameters subject to perturbations. These results 
have also been tailored to solve the problem when the state space representation of the system rather 
than transfer function is given. In this case, some entries of state space matrices are considered as 
perturbed parameters. 

Moreover, the techniques developed here overcome the difficulties that are encountered by 
the traditional frequency domain classical domain analysis and design when multiple parameter 
variations are presented. We have successfully developed new tools that enable us to construct Bode, 
Nyquist, and Nicholes templates for the control system containing multiple parameter variations. 
Using these innovative techniques, designers now have access to information which is vital to the 
design of control systems. As a result, the Matlab based toolbox for Robust Parametric Control 
has been implemented, and the final stage of testing is being conducted. A portion of this toolbox 
is available through a Mathworks 5 anonymous ftp site for those interested in this research. The 
results have also been extended to the problem of system identification. Our technique effectively 
constructs a parameterized transfer function from the given set of input - output pairs of data that 
represents the behavior of the system to be modeled. 

For the past three years of this NASA support, numerous important publications have been 
generated. Our recent book, Robust Control: The Parametric Approach published by Prentice Hall 
Publishing Company, contains the details of this research. 

The remainder of this report is organized as follows: In section 4, we describe the basic math- 
ematics used in robust stability problems for interval systems. These results are further extended 
to achieve frequency domain properties of such systems. In conjunction with this frequency domain 
information, design and analysis techniques which are similar to those used in classical control are 
described in detail. In section 5, the result of Section 4 is improved to deal with more complex sys- 
tems (multilinear interval systems). All the results developed in the framework of interval systems 
are extended to this class of systems. The main tool we have employed here is the Mapping Theo- 
rem. Section 6 shows an interval system modeling technique, and it is successfully applied to a large 
space structure with variable load. This section demonstrates successful application of the results. 
Finally, using the results obtained for multilinear interval systems, some techniques are presented 
which determine robust stability of the system when the state space representation is given. 
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4 ROBUST CONTROL VIA INTERVAL SYSTEMS 


Frequency response methods play a fundamental role in the fields of control, communications and 
signal processing. Classical control focuses on the frequency domain properties of control systems and 
has developed design methods based on simple but powerful graphical tools such as the Nyquist plot, 
Bode plots, and Nichols Chart. These techniques are well known and are popular with practicing 
engineers. However, they were developed for a fixed nominal system and in general are inapplicable 
when several uncertain parameters are present. In these situations it is necessary to evaluate the 
frequency domain behaviour of the entire family of systems in order to effectively carry out analysis 
and design. 

A brute force approach to this problem (grid the uncertainty set) can be avoided by assuming 
a certain amount of structure for the perturbations even if such an assumption introduces some 
conservatism. In this section we shall consider the class of linear interval systems where the uncertain 
parameters lie in intervals and appear linearly in the numerator and denominator coefficients of the 
transfer functions. For example, the family of transfer functions 

4s 3 -I- a 2 s 2 + ai s + 5 
[S) ~ s 4 T10/335 3 + /325 2 + (A+2 7i )s 

where a 2 , a i,/? 3 , /? 2 ,/?i,7i vary in independent intervals is a linear interval system containing six 
interval parameters. In this example, the uncertainty template G (ju) at each frequency u) is a 
complex plane set generated by the parameter vector ranging over the six dimensional parameter 
box. With the results to be developed in this section we will be able to replace G(s) by a subset of 
systems Ge(s). This extremal subset will allow us to constructively generate the exact boundary of 
the uncertainty template by means of a set of one parameter problems. These extremal systems will 
allow us to exactly calculate the boundaries of the Bode, Nyquist and Nichols plots of all transfer 
functions in the control system. They also can be used to calculate the worst case gain, phase, and 
parametric stability margins over the uncertain set of parameters. The utility of these concepts in 
control system design is illustrated by giving examples which robustify classical design techniques 
by incorporating parametric uncertainty. 

We begin by considering an interval plant connected to a fixed feedback controller and develop 
the appropriate mathematical machinery for this system. The generalized Kharitonov segments serve 
to define the extremal systems. Using these systems, we calculate the boundaries of the image sets of 
various system transfer functions evaluated at s = ju>. These include the characteristic polynomial, 
open and closed loop transfer functions, sensitivity and complementary sensitivity and disturbance 
transfer functions. We also evaluate the worst case stability margins using these extremal systems. 
These results depend on some simple geometric facts regarding the sum and quotients of complex 
plane sets. We then generalize these results to the larger class of linear interval systems using 
essentially the same geometric ideas. 

4.1 INTERVAL CONTROL SYSTEMS 

Consider the feedback system shown in Figure 1 with 



Figure 1: A unity feedback interval control system 


F(s) := 


Fi(s) 

F2(s)' 


G(s) := 


N(s) 

D(s)- 


( 1 ) 


We suppose that F(s) is fixed but G(s) contains uncertain real parameters which appear as the 
coefficients of N(s) and D(s). Write 


D(s) := CLq + 0>i S + Cl2 s ^ + G 3 S 3 + • * • + a n _is n 1 + 

N (s) bo + T &2 fi2 T &3S 3 * • • + b m — \s m 1 + b m s m (2) 
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where a* E [a k ,a£], for k E n {1, • • • ,n} and 6* € [5* for k E m. Let us define the interval 
polynomial sets 


D(s) := {D(s) : a 0 + a\S + a 2 s 2 -{ b a n s n , a* E [a fc , a£], for ken] 

N(s):={N(s) : 6 0 + M + &2« 2 d b b m s m , b k E [&*, 6*], for fc E m} 

and the corresponding set of interval systems : 


G(s) := { : € (n(*)xd(*)) } • 


(3) 


We refer to the unity feedback system in Figure 1 as an interval control system. For simplicity, we 
will use the notational convention 

G <*> ■ 5 $ (4 » 

to denote the family (3). The characteristic polynomial of the system is 


6(s):=F l (s)N(s)+F 2 (s)D(s) (5) 

and the set of system characteristic polynomials can be written as 

A(s) F l (s)N(s) + F 2 (s)D(s). (6) 


The control system is robustly stable if each polynomial in A(s) is of the same degree and is Hurwitz. 

Let /C n (s) and ICd(s) denote Kharitonov polynomials associated with N(s) and D(s), and let 
Sjv(s) and Sd(s) denote the corresponding sets of Kharitonov segments. Recall that these segments 
are pairwise convex combinations of Kharitonov polynomials sharing a common even or odd part. 
Define the extremal subsets , using the above notational convention: 


Ge(s) := irj-l u (extremal systems) (7) 

Sd{8) K,d{s) 

Gk(s) := (Kharitonov systems). (8) 

ICd{s) 


We shall say that F(s) satisfies the vertex condition if the polynomials Fi(s) are of the form 


Fi(s) := (a iS + bi)Ui(8)Ri(s) y i = 1,2 (9) 

where U are nonnegative integers, a,, b{ are arbitrary real numbers, Ui(s) is an anti-Hurwitz poly- 
nomial, and Ri(s) is an even or odd polynomial. We recall the result given by GKT. 

Theorem 1. The control system of Figure 1 is robustly stable that is stable for all G(s) € G (s) if 
and only if it is stable for all G(s ) € Ge(s). If in addition F(s) satisfies the vertex condition , robust 
stability holds if the system is stable for each G(s) E Gk(s). 

The GKT thus reduces the problem of verifying robust stability over the multiparameter set 
G(s) to a set of one parameter stability problems over Gjj(s) in general, and under the special 
conditions on F(s) stated, to the vertex set G/c(s). In the rest of this section we shall show that 
the systems Ge(s) and Gk(s) enjoy many other useful boundary and extremal properties. They 
can be constructively used to carry out frequency response calculations in control system analysis 
and design. In fact, it will turn out that most of the important system properties such as worst 
case stability and performance margins over the set of uncertain parameters can be determined by 
replacing G(s) E G(s) by the elements of G{s) E Ge(s). In some special cases one may even replace 
G{s) by the elements of Gk (s). The results are first developed for interval plants for the sake of 
simplicity. They hold for the more general class of linear interval systems as indicated in section 4.5. 

4.2 FREQUENCY DOMAIN PROPERTIES 

In order to carry out frequency response analysis and design incorporating robustness with respect 
to parameter uncertainty we need to be able to determine the complex plane images of various 
parametrized sets. In this section we will develop some computationally efficient procedures to 
generate such sets. We shall first consider the complex plane images of A (s) and G(s) at s = ju. 
These sets, called uncertainty templates, are denoted A (jw) and G(ju). Since N (s) and D(s) are 
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interval families, N(ju;) and D(jcj) are axis parallel rectangles in the complex plane. F\ (ju)N(ju) 
and F 2 {ju)D{ju) are likewise rotated rectangles in the complex plane. Thus A (ju) is the complex 
plane sum of two rectangles whereas G(ju;) is the quotient of two rectangles. We assume here that 
0 £ D(ju>). If this assumption fails to hold we can always “indent” the ju axis to exclude those 
values of u which violate the assumption. Therefore, throughout this section we will make the 
standing assumption that the denominator of any quotients exclude zero. 

The next lemma will show us how to evaluate the sum and quotient of two complex plane 
polygons Qi and Q 2 with vertex sets Vi and V 2 i and edge sets E\ and E 2 , respectively. Let <9( ) 
denote the boundary of the complex plane set (•). 

Lemma 1. (a) d(Q 1 + Q 2 ) C (Vi + E 2 ) U (E 1 + V 2 ), (b) d (^) C f U 

This lemma shows us the interesting fact that the boundaries of sums and quotients of two 
polygons can be determined by sums and quotients of the corresponding vertex-edge pairs. Figure 2 
illustrates that the boundaries of the sum of two four sided polygons are obtained by generating the 
sum of all segment-edge pairs. Similarly, Figure 3 shows the sum of two general polygons. 

It can be shown that the inverse image of a line segment which excludes the origin, is an 
arc of a circle passing through the origin. This is shown in Figure 4. Therefore the inverse image 
of a polygon is bounded by such arcs as shown in Figure 5. To determine A (ju) and G(ju) we 
note that the vertices of N (ju;) and D(ju>) correspond to the Kharitonov polynomials whereas the 
edges correspond to the Kharitonov segments. The set of points K,s(ju) are therefore the vertices 
of N (ju) and the four lines Sw(ju) are the edges of N(ju;). Fi(ju)N(ju) is also a polygon with 
vertices Fi(ju)Kx(ju) and edges F\(ju)SN(ju). Similarly, F 2 (ju)/Co(ju) and F 2 (ju)So(ju) are 
the vertices and edges of the polygon F 2 (ju)D(ju ;). The ju image of the extremal systems G E (s) 
defined earlier exactly coincides with these vertex-edge pairs. Let 

(N(s)xD(s)) E := (fC N (s)xS D (s)) U (S N (s)x/C D (s)). (10) 

Recall that the extremal systems are 

GeW: ={^ : (mCW)€(NMx D W ) B } ; =|^u|^ (11) 

and define 

A e (s) := {Fi(s)N(s) + F 2 (s)D(s) : (N(s),D(s)) G (N(s)xD(s)) E } . (12) 

We cam now state an important result regarding the boundary of image sets. 

Theorem 2. (Boundary Generating Property) 

a) dA (ju) C A E (ju) b) dG (ju) C Ge (ju) 

Example 1. Consider the problem of determining the frequency template of the interval plant 

hs + bo 

Cj(s) — t/ . — ^ 

d(s) a 2 S l +a\s + ao 

where the parameters vary as follows: 

ao € [1,2], a\ G [2,3], a 2 G [2,3], 6o € [1,2], b\ G [2,3]. 

The Kharitonov polynomials of d(s) and n(s) are: 

K^(s) = 3s 2 + 2s + 1, Kj(s) = 3s 2 + 3s + 1, K$(s) = 2s 2 + 2s + 2, Kj(s) = 2s 2 + 3a + 2 

and 

K\(s) = 2s + 1, Kfo) = 3s + 1, K*(s) = 2s + 2, K*(s) = 3s + 2. 

Thus, the boundary of the entire frequency domain template is obtained by the frequency evaluation 
of the following 32 systems: 

KU f) A Kjjs) + (1 - A )K*(s) 

\K j d {s) + (1 - \)K*(s) ’ K^s) 


for 


i = 1,2, 3,4; (j,k) G {(1,2), (1,3), (2, 3), (3, 4)}. 


Figure 6 shows the template G(jtj) at u = 1. 
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Figure 5: Polygon and the inverse of polygon 
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Figure 6: Frequency domain template G(ju?) (Example 1) 


4.2.1 Closed Loop Transfer Functions 


Referring now to the control system in Figure 1, we consider the following transfer functions of 
interest in analysis and design problems: 


yQ) 

u(s) 


= 0 ( 8 ), 


T°(s) := ^ = F ( S ) G ( S )’ 

e(s) 


T u (s) := 


F(s) 


u(s) 

r(s) 1 + F(s)G(s) ’ 


u(s) 

e(s) 


= F(s) 


i •= = 

1 

r(s) 

1 + F(s)G(s) 

i y( s ) _ 

F(s)G(s) 

r(s) 

1 + F(s)G(s) 


(13) 


As G(s) ranges over the uncertainty set G (s) the transfer functions T°(s), T v (s), T u (s), T e (s) range 
over corresponding uncertainty sets T°(s), T w (s) ) T u (a) , and T e (s), respectively. In other words, 


T°(s) := {F(s)G(s) 


rj>U 



F(s) 

1 + F(s)G(s) 


G(s) G G(s)}, T*(s) := | r — : G(s) G G(s)j (14) 

G <*> 6 G( *>} • T "W : = : G « 6 G <»>} ■ 


We will now show that the boundary generating property of the extremal subsets shown in Theorem 2 
carries over to each of the system transfer functions listed above. In fact, we will show that the 
boundary of the image set at s — , the Nyquist plot and Bode plot boundaries of each of the above 

sets are all generated by the subset Ge (s). Introduce the subsets of (14) obtained by replacing G(s) 
by G e (s): 

T|(s) := {F(s)G(s) : G(s) G G E (s)}, T|(s) := { l + F ^ )G{s) ■ G(s) G G E (s)}(15) 

T '<*» : = { 1 + : G(>> 5 GeW } ' W ; = { 1 + F(°G(a) ; G <*> e G ' W } • 

The main result can now be stated. 


Theorem 3. For every uj > 0, 

(a)dT°(jLj) C T |0'w), (b)dT'(juj) C T| (Jcj), (c)dT u (ju>) C TS(jo>), (d)dT»(jw) C T v E (ju>) 
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Remark 1. This result shows that at every u) > 0, the boundary of the image set of each transfer 
function in (14) is contained in the corresponding image set of the extremal systems. We point out 
that in the definition of the interval plant G(s) we assumed that the numerator and denominator 
parameters were independent. This assumption obviously does not hold any longer when we deal 
with say, T y (s) where the numerator and denominator depend on some common parameters. It 
is therefore useful to know that the boundary generating property of the set Ge (s) carries over 
nevertheless. In a later section we will show that the boundary generating property of the set Ge(s) 
will hold much more generally. 


Example 2. Consider the system given in Example 1. Let us assume that F(s) = 1 and that we 
wish to calculate T y (juj). 


T ' < ’ )= {irGw :GWeGW }' 


For this example, we have 


T*(s) = 


b\ s + bo 

d 2 g2 4 (&i 4 b\)s 4- (ao 4- bo) 


where a 0 € [1,2], a x 6 [2,3], a 2 6 [2,3], b 0 € [1,2], b x € [2,3]. 

The denominator and numerator polynomials are dependent on some of the same perturbing 
parameters. However, Theorem 3 handles this dependency. Since dT y (jui) C Tg(ju;), it is enough 
to construct the following template: 




In other words, we simply replace G(s) by Ge(s) in T y . Thus we need to construct the image of 
the following transfer functions for A € [0, 1]: 

KL ju) \K*(ju) + (1 — X)K l n (jio) 

A KS(jcj) + (1 - A )K l d (juj) + K'Jju) ’ A K^ju) + (1 - A 4 K'Jju) 

for i = 1,2, 3, 4 and (kj) € {(1, 2), (1, 3), (2, 3), (3, 4)}. The sets of Kharitonov polynomials corre- 
sponding to the denominator and numerator polynomials are defined in Example 1. The frequency 
domain template at cj = 1 is shown in Figure 7. 



Figure 7: Frequency domain template T v (ju) (Example 2) 


11 




4.3 NYQUIST, BODE, AND NICHOLS ENVELOPES 

In the previous section we established that at a fixed frequency the image set template of various 
transfer function set can be constructed from the corresponding extremal transfer function set Ge(s). 
In control system design it is important to see the behaviour of these transfer function families for 
all frequency. Denote the Nyquist plot of the family G (s) as 


G = Uo<a;<ooG(jU^). (16) 

The boundary of G is the Nyquist envelope of G. Similarly, the Nyquist plots of T °(s), T e (s), 
T u (s), and T y (s) are denoted respectively by T°, T e , T u , and T y . Prom the boundary property 
of Theorem 3, it follows that the envelopes of these Nyquist plots are generated by the extremal 
systems. 

Theorem 4. (Nyquist Envelope) 

The Nyquist plots of each of the transfer function sets T°(s), Tt'(s), T u (s), and T e (s) are bounded 
by their corresponding extremal subsets: 

(a)d GcG e , (b)d T°CT° E) (c)dT»cT y E , (d)d T u C Tg, (e)OT e C T|. (17) 


The Bode envelopes of each of the system transfer functions G(s), T°($), T e (s), T u (s), and T v (s) 
can also be constructed if we can determine the maximum and minimum values of the magnitude and 
phase of the family of systems at each frequency. From the boundary relations given in Theorem 3, 
it again follows that these maximum and minimum values occur over the subset Ge(s). This leads 
us to the following result. 

Theorem 5. (Bode Envelope) 

The Bode magnitude and phase envelopes of each of the transfer function sets G(s), T°(s), T e (s), 
T u (s), T y (s) are generated respectively by the extremal subsets, Ge(s), T|(s), T|(s), Te(s), Tg(s). 

It is instructive to interpret this result from the geometric point of view. Consider the Bode plot of 
the interval family G(s). For a fixed to * , the image of the interval transfer function is the quotient of 
two convex polygons in the complex plane where the polygons represent the images of the numerator 
and denominator polynomial families. This is depicted in Figure 8. 



Figure 8: Extremal magnitude and phase of N and D (ju*) 
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From Figure 8(a) and (b), we can see that the maximum and minimum magnitudes of G(ju*) 
occur on one of the vertex-segment combinations of these polygons. Since the parameters in the 
numerator and denominator are independent we have: 


(a) 


max 




max |N(jo;*)| N(ju>*) 

min|D(jw*)l D(jW) 


min |N(ju;*)| 
max |D(jcj*)| 


While the maximum magnitude always occurs at a vertex the minimum can occur on an edge. The 
maximum and minimum points will generate the extreme points of the Bode magnitude envelope at 
the frequency u*. On the other hand, as shown in Figure 8(c) and (d), the extreme points on the 
phase envelope are always generated by vertex- vertex pairs: 


(C) 

maxarg 

(d) 

min arg 


n (jun 

D(jw*) 

D(jw*) 


— maxargN(ju;*) — minargD(ju;*) 


= minargN (jco*) — maxargD(ju;*). 


The relations given in the Theorem above are useful in constructing the Bode magnitude and 
phase envelopes as well as Nyquist envelopes. In classical control design techniques, the Nichols 
Chart is also a popular computational tool. For a fixed transfer function, the Nichols Chart is a 
plot of the magnitude versus the phase with frequency as a running parameter. When a family 
of parametrized systems is involved, we get at each frequency a magnitude-phase template. When 
this template is swept over all frequencies, we get the Nichols Envelope of the family. By using the 
boundary property of Theorem 3 it can be seen that the Nichols Envelope is also generated by the 
extremal systems. 

Theorem 6. (Nichols Envelope) 

The Nichols Envelope of each of the transfer function sets G (s), T °(s), T e (s), T u (s), T y (s) are 
generated respectively by the extremal subsets, Ge(s), Tg(s), T|(s), Tg(s), T|(s). 

The Nichols template is obtained by mapping the points of the Nyquist plane (see Figure 9) into 
the Nichols plane (see Figure 10). This gives the exact Nichols template. 



The Nichols envelope can also be generated approximately from the Bode magnitude and 
phase envelopes of the family. At each frequency we draw a rectangle in the magnitude-phase plane 
(Nichols plane). The dashed rectangle shown in Figure 10 corresponds to the magnitude and phase 
ranges obtained from the Bode envelopes. When these overbounded images of G (juj) templates 
are swept over frequency, we obtain the approximate Nichols envelope which contains the actual 
envelope. We illustrate these computations with an example. 


13 



90°= e 2 



Example 3. Consider the plant and controller 


G(s) = 


b\8 + bo 
a 2 s 2 + ais + ao 


C(s) = 


8 2 + 2s + 1 
a 4 + 2s 3 + 2s 2 + a 


where the plant parameters vary as 

6i € [0.1, 0.2], 6 0 e [0.9, 1.1], a 2 € [0.9, 1.0], ai G [1.8, 2.0], oo € [1.9, 2.1]. 
Figures 11, 12, and 13 show the frequency domain plots for this example. 
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4.3.1 Conservatism of the Envelopes 

For a fixed system, all of the above three frequency plots provide the same information on the 
system. However, for the case of a parametrized family of systems, the situation is quite different. 
The Nyquist and Nichols plots of a fixed system can be regarded as a string in the Nyquist or Nichols 
plane. For a parametrized family of systems, the Nyquist or Nichols plots therefore consists of a 
family of strings. However, the envelope of the plot is, in general, not a string that belongs to the 
family. In other words, there is no system in the family which generates the entire boundary of the 
envelope itself. On the other hand, every point on the boundary of the envelope has a string passing 
through it as shown in Figures 14 and 15. 



Figure 14: Nyquist envelope and Nyquist plots of individual member systems 



Figure 15: Nichols envelope and Nichols plots of individual member systems 

The Bode envelopes in fact correspond to an overbound of the image of the complex plane set 
G(ju;) at each frequency u>. Considering Figure 16, we also see that each point on the boundaries of 
the envelopes came from a true member system of the family. From the Bode envelopes, we notice 
that the smallest gain margin of the family is K. However, the member system passing through 
the point in the magnitude plot does not correspond to the phase crossover point in the 
phase plot. Due to this phenomena, the true smallest gain margin might be bigger than K. A 
similar argument can be made for the case of the smallest phase margin (see Figure 17). This 
is because of the independent evaluation of the magnitude and phase envelopes of the family. In 
other words, despite the fact that each point on the boundaries of the Bode magnitude and phase 
envelopes, comes from some parameter in the family, the latter envelopes taken jointly, represent 
only the approximate set which corresponds to the dashed box in Figure 10 and equivalently, the 
dashed portion of the disc in Figure 9. Therefore, the smallest gain and phase margins read from 
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the Bode envelopes would be conservative, with the degree of conservatism depending upon how big 
the actual image inside the dashed box is (see Figure 10). 

Remark 2. The above boundary results remain valid if Fi(s) are complex functions rather than real 
polynomials. This is useful in applications. For example, in systems containing time-delay Fj(s) 
could contain terms such as e ~ T * . Also many robust performance problems reduce to verifying 
robust stability under real parameter perturbations using complex compensators. 

4.4 EXTREMAL STABILITY MARGINS 

In this section, we deal with the calculation of guaranteed stability margins. We will first consider 
the gain and phase margins at the loop breaking point “m” shown in Figure 18. When G(s) 


TO 


m 


F(s) 


G (s) 


Figure 18: A unity feedback system 

contains parameter uncertainty the worst case values of gain and phase margins over the parameter 
uncertainty set axe important measures of robust performance. We show that these worst case 
margins in fact occur over the subset of extremal plants Ge(s). 


4-4-1 Guaranteed Gain and Phase Margins 

Suppose that F(s) robustly stabilizes the family G(s). The gain margin of the system for a fixed 
G(s) € G(s) is defined to be the smallest value Kq for which 

(1 + K)F l (8)N(s)+F 2 (s)D(8) 
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remains stable for all K € [0, Kq)- Similarly, the phase margin of the system for a fixed G(s) € G(s) 
is defined to be the smallest value 6g for which 


e j0 F 1 (a)N( 8 ) + F 2 (s)D(s) 


remains stable for all 9 € [0, 9q)- The worst case gain and phase margins are: 


K * := inf K G 

G(»)e G(») 


6* inf 0 g - 

G(»)e G(») 


(18) 


Theorem 7. (Extremal Gain and Phase Margin) 

I) 

K* = inf K g , 9* 

G(«)GGe(«) 


inf 8g • 

G(s)€Ge(») 


II) If Fi(s) are real and satisfy the vertex conditions specified in part II of GKT , then we have 

K* = inf K g 

G(8)€G k (8) 


The proof of this theorem readily follows from the fact that A(jcj) and AeC?^) share identical 
boundaries. Moreover, when Fi(s) satisfy the vertex conditions the proof follows from the fact that 
AW is stable if and only if Ak(s) is stable. 


4-4-2 Worst Case Parametric Stability Margin 

We now consider the worst case parametric stability margin. We assume as before that F(s) robustly 
stabilizes G(s). The parameter p of dimension t consists of the coefficients of N(s) and D(s) and 
varies in the hypercube IT. We will rewrite S(s) in (5) as S(s, p) to emphasize its dependence on 
the parameter p. Let n E denote the parameter subset of II corresponding to the extremal systems 
Ge(s). Let n K denote the parameter subset of II corresponding to the extremal systems G K (s)- 
Let || • || denote any norm in R* and let V u denote the set of points u in R ( for which S(s, u) 
is unstable or loses degree (relative to its generic degree over II). Let 

p( p) = inf Up - u llp 

uev u 

denote the radius of the stability ball (measured in the norm || • ||) and centered at the point p. This 
number serves as the stability margin associated with the point p and indicates its distance from 
instability. If the box II is stable we can associate a stability margin, denoted p(p), with each point 
in II. A natural question to ask then is: What is the worst case value of the parametric stability 
margin over II and what is the point where it occurs? An answer to this question gives an indication 
of how close one can get to instability over the box II. 

Define a mapping from II to the set of all positive real numbers: 

n -A ft+\{ 0 }, p— >p(p) 

Our question stated in terms of functions is: Has the function p(p) a minimum and is there a precise 
point in II where it is reached? The answer is provided by the following theorem: 

Theorem 8. (Extremal Parametric Stability Margin) 

I) The minimum value over II of p(p) is reached at a point on the extremal set IIe . 

II) If Fi(s) satisfy the vertex conditions, the minimum value of p(p) is reached at a point on the 

extremal set IIk- 


Example 4- 


Consider the unity feedback control system with 
ais + u-o 


G{s) = 


b 2 s 2 4- b\S + 6q 


and F(s) = 


s 2 4- 2s 4- 1 
+ 2s 3 4* 2s 2 4- s 


where 


a o € [1 — €, 1 + e], 

b\ £ [1.9 — e, 1.9 4- c], bo € [2 — c, 2 4* c]. 


a\ £ [0.15 -e, 0.15 + c], 
b 2 E [0.95 — e, 0.95 e], 
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We wish to find the largest excursion of parameters allowed or equivalently the maximum value of 
e for which closed loop stability is preserved. The image set plot at c = 0.128, shown in Figure 19, 
reveals that the phase difference reaches 180° at e w 0.128 (Figure 20). It is clear that any value 
of e smaller than 0.128 results in smaller image sets than the ones shown in Figure 19. Moreover, 
if we sweep the frequency the connected image sets will form an envelope which does not contain 
the origin for e = 0.128. Therefore, an envelope corresponding to any value of e smaller than 0.128 
cannot contain the origin. Thus, we conclude that the feedback system remains stable under the 
perturbations bounded by this maximum value of c. 



Real 


Figure 19: Image set plot at e = 0.128 (Example 4) 



Figure 20: Phase plot (Example 4) 


4.5 LINEAR INTERVAL CONTROL SYSTEMS 

The results given so far in this section assume that G(s) is an interval plant, namely the ratio of 
interval polynomials. Each of the results given shows that a particular calculation involving G() 
can be replaced by the corresponding calculation over the one parameter subsets Ge()* This sort of 
simplification actually carries over to the more general class of linear interval systems G (s) defined 
below. Instead of repeating all the previous results, we show in this section how to construct the 
extremal subsets Ge(s) and leave the details of the proofs to the reader. 

We still consider the unity feedback system shown in Figure 1 with 


F(s) := 


Fi(s) 

F2(sV 


G(s) := 


**{») 

D(*y 


(19) 
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We suppose that F(s) is fixed, <?($) is subject to parametric uncertainty and is now of the form 
N(s) := Li(s)A\(s) 4- L 2 (s)A 2 (s) + • • • + L u {s)A u {s) 

D(s) := M\(s)Bi {a) + M 2 (s)B 2 (s) + • • ■ 4 M v (s)B v (s) (20) 

with i4i(«) and Bj(s) being independent interval polynomials. Write 

Ai(s) := a l 0 4 a\s 4 a\s 2 4- a l 3 s 3 4 (- 4- a ni s ni 

Bj(s) := 4 + b{s + V 2 s 2 + bis 3 + --- + (21) 

where a\ G [a). _ ,a). + ], for k 6 n* and b 3 k G [6j - ,6j + ], for k G dj, and Lj(s) and M,(a) are fixed 
polynomials in s . Let 

AW := [A! W, AaW, ' • • , A«W], flW == [*i W, ftW. • • • , *.W]- (22) 

Let us define the sets 

Aj(s) := {Aj(s) : a' 0 + ajs + a 2 s 2 4 b a},.s ni , a\ € [ajT,a* + ]> f° r * e Hi) 

Bj(s) := {^(s) : ftj + frjs + i^s 2 + • • • + fcj.A fej G (4 >*44 for /c G dj } 

N W : = j^^iW^iW : (^i( s )>-'->^u(s)) € Ai(s) x ••• x A u (s) 

U=1 

B(s) := lj2Mi(s)Bi(s) : (Bi(s), ■ ■ ■ ,B v (s)) G Bx(s) x • • • x B„(«) 
l*=l 

and the corresponding set of uncertain systems 

G W := {^J : (N(s),D(s)) G (N(s)xD(s))} := ® (23) 

is called a linear interval system. The Kharitonov polynomials and segments associated with A;(s) 
are denoted /C^(s) and respectively. K, Bh {s) and SB h (s) are defined similarly. Now let us 

define 

Ak(s) := ICai {s) x • * * x ICa v («), B K (s) := Kb x (5) x • * • x ICb v (5) (24) 

and 

AeW := KaM X X £a*_xW X S Ak {s ) X )C Ak+1 (s)--- x £,t u (s), 

Bg(s) := IC Bl (s) x • • • x £ b> _, (s) x S Bj (s) x JC Bj+1 W • • • x /C B „ (s) (25) 

for k £ u and j 6 v. Let 

AeW := UttAfeW, B e (s) := U^B^s). (26) 

Now introduce 

NkW := {n(s) : N(s) = L i (s)A i (s), (A t (s), • • • , A u (a)) G A K (s)| 

N e (*) := {N(s) : N(s) = £ Li(s)Ai (*), (A^e), ■ ■ ■ , A u (s)) G A E (s)} • (27) 

Similar definitions hold for Dk(s) and De(s), and the extremal subsets are: 

(N(s)xD(s)) E := (Nk(s)xD e (s)) U (N E (s)xD K (s)). (28) 

The extremal subset of the family G (s) is: 

: = {£$ ; <"<*>' D <*» 6 == < 29 > 
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In words the extremal set is obtained by fixing all the A { , Bj at Kharitonov vertices except one 
and letting this one range over the Kharitonov segments. We remark here that in the more general 
case where some of the A{ and Bj happen to be identical the same procedure for constructing the 
extremal subset will work except that this constraint must be imposed on Ge(s). 

With the above definitions we can verify that N (ju) and D(jcj) are polygons. The vertices of 
’N(ju) are contained in the set Nk ( jcj) and the edges of N (ju) are contained in Ne (ju)- Similar 
relations hold for D (juj). From this, it follows easily that for the class of linear interval systems 
G(s) defined in (20) - (23) each of the statements given in Theorems 1-8 remains valid with 
Ge(s) defined as in (29). The conditions are constructive because here again Ge(s) is a set of one 
parameter families. We illustrate the above calculations in the example below. 

Example 5. Consider the linear interval system 

G( , _ Li(a)i 4 i(a) 5(7 i a + 70) 

M 1 (s)Bi(s) + M 2 (s)B 2 (s) s(a 3 s 3 + a 2 s 2 ) + (/3 2 s 2 + fas) 

where 71 € [0.40,1.60], 7o € [19.40,20.60], a 3 € [0.40,1.60], a 2 G [7.40,8.60], p 2 G [31.4,32.6], 
/?! E [74.45,75.65]. Here, we illustrate the construction of Bode, Nyquist, and Nichols templates. 
First we write the Kharitonov polynomials of Ai(s), 2?i(s), and ^(s): 

K\ x {s) = 1.6s 3 + 8.6s 2 K 2 m (s) = 0.4s 3 + 8.6s 2 K 3 Al {s) = 1.6s 3 + 7.4s 2 

K* Ai (s) = 0.4s 3 + 7.4s 2 ATgj (s) = 32.6s 2 + 74.45s K 2 Bi (s) = 32.6s 2 + 75.65 
K 3 Bl (s) = 31.4s 2 + 74.45s K% t (s) = 31.4s 2 + 75.65s K 1 ^ (s) = 0.4s + 19.4s 

K% 2 (s) = 1.6s 2 + 19.4s K 3 B2 (s) = 0.4s 2 + 20.6s K^js) = 1.6s 2 + 20.6s 

Then the extremal systems we will deal with are 

Inis) (AKj ti (s) + (l-A)j<:* i (s)) 

M 1 {s)K l Bl ( S ) + M 2 (s)K 1 B2 (s) ’ 


Q LiisJK^js) 

M,(s) (A K k Bi (s) + (1 - X)K l Bl (s)) + M 2 (s)K j B2 (s)’ M 1 (s)Ki i (s) + M 2 (s) (A K^(s) + (1 - A )K l Ba (s)) ’ 

for i,j = 1,2,3, 4; ( k,l ) G {(1, 2), (1, 3), (2, 4), (3, 4)}. Figures 21, 23 and 22 show the Nyquist, 

Bode and Nichols envelopes of the system respectively. 



Figure 21: Nyquist templates of linear interval system (Example 5) 
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Figure 22: Nichols templates of linear interval system (Example 5) 




rad/sec 

Figure 23: Bode envelopes of linear interval system (Example 5) 
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4.6 POLYTOPIC SYSTEMS 

In many models containing parametric uncertainty the only convenient representation is one where 
each independent parameter is explicitly displayed. In this case, we write 


aiLi(s) T 02 ^ 2 ( 5 ) + • • * + a u L u (s) 
W b 1 M 1 (s) + b 2 M 2 (s) + --- + b v M v {s) 

(30) 

and let 

G(a) := (G(a) : a," < a t < aj, bj < bj < bj, i £ u, j € v] 

(31) 

denote the family of uncertain systems. Let 


P * — [&1 ? ^2 j * ' ‘ ) &u 1 ^ 1 1 ^2 > * * * 1 

(32) 

denote the parameter vector and 


II := {p : aj < o* < aj, bj < bj < i € u, j € v} 

(33) 

the uncertainty polytope. Write (30) a s G(s, p) to emphasize its dependence on p. For convenience 
we refer to the family of systems G(s, p) : p G II represented in the form (30) as a polytopic system. 
Of course a polytopic system is just a special case of a linear interval system with the uncertain 
polynomials being of degree zero. 

Let IIe denote the exposed edges of II and introduce the extremal systems 

G E (a) := {G(a,p) : pen E }. 

(34) 


Now suppose that s* is an arbitrary point in the complex plane, where the image set G(s*) needs 
to be found. In the Hurwitz case we might have s* = ju and in a discrete time control system, for 
example, s* could be a point on the unit circle. We have the following result. 


Theorem 9. 


dG(s*) C G e (0. 

All the results of Theorems 1 to 8 carry over to this general case with the extremal set Ge(s) being 
defined by (34). The more general case where some a*, bj are identical can be handled by imposing 
the same constraint on Ge(s). 


4.7 LINEAR FRACTIONAL TRANSFORMATIONS OF INTERVAL SYSTEMS 


We have thus fax assumed that our uncertain system is described by an interval system, a linear 
interval system, or a poly topic system model. In each of these cases we developed extremal and 
boundary results under the assumption that the numerator and denominator parameter sets are 
disjoint or that they perturb independently. On the other hand, we saw that these extremal and 
boundary results carry over to all the feedback system transfer functions even though, in some 
of these, the numerator and denominator do contain common uncertain parameters. A natural 
question to ask therefore is: What is a useful general model with interval parametric uncertainty, 
incorporating dependencies between numerator and denominator parameters, for which boundary 
or extremal results can be obtained? It turns out that a very broad class of systems of this type can 
be encompassed under the class of linear fractional transformations of an linear interval system. 

Therefore, let us suppose that G(s) is a linear interval system (see (23)) with independent 
uncertain interval polynomials and let Ge(s) denote the extremal subset as defined in 

Section 4.5. Let P(s), Q(s), R(s), S(s) be four arbitrary fixed functions and consider the transfor- 
mation 

P(a)G(a) + Q(a) 
w R(s)G(s) + 5(a) ‘ 


We suppose that the transformation is well defined ( R(s)G(s ) + S(s) ^ 0) and refer to such a 
transformation as a linear fractional transformation. Obviously H(s) contains, in general, common 
interval parameters in the numerator and denominator even though G(s) has only independent 
parameters in the numerator and denominator. As before let 


HW 


-{ 


P(a)G(a)+Q(a) 
R(s)G(s) + 5(a) 


: G(a) 6 G(a) 
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and let 


r fwgw + ow a 

HeW = \ JJMGM + S(») ' G(s) € GeW J ' 


We show below that under mild conditions on the set P(s ), Q(s), R(s) y S(s) the boundary of the set 
H(ju;) is generated by the extremal set GeO'w). 


Theorem 10. 

we have 


If 


P(ju)S(juj) - Q(ju)R(jw) 0 
dH (ju>) C H e (j'w). 


Remark 3. With this result in hand it is easy to see that all the extremal and boundary results stated 
in this section carry over to linear fractional transformations of linear interval systems. The practical 
implication of this result is that by adjusting the set P(s),Q(s), R(s) y S(s) we can account for a large 
class of interdependencies between the perturbations. It is easy to show that all the transfer functions 
that occur in a closed loop system can be shown to be linear fractional transformations of G(s) which 
satisfy the restriction stated in the theorem above. Thus it is not surprising that the extremal results 
stated in Theorem 3 hold. Finally, we point out that the boundary generating property of Ge(s) 
continues to hold for any further linear fractional transformations applied to H(s). 

Example 6. Let us consider the block diagram shown in Figure 24. 



Figure 24: A feedback system 


The closed loop transfer function of the system is 

= P 1 (a)P A {a) + P^Gja) 

u(s) 1 + [Pi(s)P z (a) + P 2 (s)} P A (a) + [Pi(s)P 3 (a) + P 2 (s)]G(s) K h 

If G(s) belongs to an interval transfer function family G($), then the family of output transfer 
functions is 

v / P!(a)P A (a) + P 1 (a)G(a) 

W \l + lPi(s)P 3 (s) + P 2 (s)}P 4 (s) + [Pi(s)P 3 (s) + P 2 (s)]G(s) ■ K) 

P(.) Q(i) 

+ . rf ) 

[Pi(a)P 3 (a) + P 2 (s))G(s) + 1 + [Pi(s)P 3 (s) + P 2 (s)]P 4 (s) ’ V ’ 

' V ' ' v ' 

«(*) S(B) 



€ G(s)| 
€G(s)|. 


Thus, if 
we have 


P(juj)S(ju) - Q(ju>)R(ju) ± 0, 
dH(ju) C HeO'w). 


The point is, that even though H(5) is not an interval transfer function, the boundary of H(ja>) 
is captured by replacing G(5j by the elements of the extremal set G #(s) which is a considerable 
saving. Similar results hold for the Nyquist and Bode envelopes, and in fact all extremal properties 
of H (s) will occur on the subset H e(*)- 
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4.8 ROBUST PARAMETRIC CLASSICAL DESIGN 


We illustrate the utility of the above tools in extending the design techniques of classical control 
to systems containing parameter uncertainty. The requirement of robust design is that the design 
specifications must be satisfied over the entire parameter set. Thus the worst case values must be 
acceptable. Since these worst case values occur over the extremal set Ge(s)> it suffices to verify that 
the specifications are met over this set. 


4.8.1 Guaranteed Classical Design 
Example 7. Consider the interval plant 

c( v ^0 

d 3 s 3 + d 2 s 2 + d\s + do 


with n 0 € [10,20], d 3 G [0.06,0.09], d 2 £ [0.2, 0.8], d x G [0.5, 1.5], d 0 = 0. The objective is to design 
a controller so that the closed loop system 

a) is robustly stable under all parameter perturbations, 

b) possesses a guaranteed phase margin of at least 45°. In other words the worst case phase margin 

over the set of uncertain parameters must be better than 45°, and 

c) possesses a bandwidth greater than or equal to 0.1 [rad/sec] with a reasonable value of resonant 

peak M p . 

We simply follow the standard classical control design techniques with the new tools developed here. 
First we construct the Bode envelopes of the open loop interval plant G(s). From the magnitude 
and phase envelopes, we observe that the worst case phase margin of 50° is achieved if the gain 
crossover frequency uj' c is at 0.5 [rad/sec]. In order to bring the magnitude envelope down to 0 [dB] 
at the new crossover frequency u ? c , the phase lag compensator of the form 


C(s) 


1 + aTs 
1 + Ts' 


a < 1 


must provide the amount of attenuation equal to the minimum value of the magnitude envelope at 

max |G(jw c )| = -201og 10 a[d£]. 

G(]<j' c )€ GOO 

Thus, we have a — 0.01. Now we choose the corner frequency to be one decade below u' c : 


1 

aT 


10 


a;' =0.5 


and we have T = 2000. Finally, the resulting phase lag compensator becomes 


C(s) = 


1 4- 20s 
1 + 2000s* 


Here we give some analysis to check whether the controller C(s) satisfies the given design require- 
ments. The robust stability of the closed loop system with the controller C(s) may be easily de- 
termined by applying Theorem 1. Moreover by GKT since we use a first order controller it is only 
necessary to check the stability of the following vertex polynomials: 


where 


(1 + 2000s) K' d (s) + (1 + 20s)A£(s), i = 1,2, 3, 4; j = 1,2 


K l n (s) := 10 
K\(s) := 0.5s + 0.8s 2 + 0.09s 3 
K%(s) := 0.5s + 0.2s 2 + 0.09s 3 


K 2 n (s) := 20 

K 2 d (s) ■- 1.5s + 0.8s 2 + 0.06s 3 
K A d {s) := 1.5s + 0.2s 2 + 0.06s 3 . 


The eight polynomials above are stable and this shows that the closed loop system remains stable 
under all parameter perturbations within the given ranges. Figure 25 shows the frequency response 
(Bode envelopes) of G(s) (uncompensated system) and C(s)G(s) (compensated system). 
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Figure 25: Bode envelopes of G(s) and <7(s)G(s) (Example 7) 


Clearly, the guaranteed phase margin requirement of 45° is satisfied. The guaranteed gain 
margin of the system is 12 db. The closed loop response |M(jcj)| called T y (jcj) in (14) is shown in 
Figure 26 where 

M(,'w) := M(a)|, =iw := {a/(s) : ^ G(s)eG( a )|. (35) 

Note that the |M(ju;)| envelope shown in this figure is calculated from the result of Theorem 4. 
Figure 26 shows that the M p of every system in the family lies between 1.08 and 1.3886. This also 
shows that the bandwidth of every system in the family lies in between 0.12 and 0.67 frad/sec]. 
Thus, the design objective is achieved. Figure 27 shows the Nyquist plot of C(s)G(s). The center 
of the M-circle in Figure 27 is 


K 

i -K 


>o I = 




and the radius of the circle is 


r = 


M p 

1 -Ml 


1.496. 
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Figure 26: Closed-loop frequency response (Example 7) 



Figure 27: Nyquist envelope of C(s)G(s) with M- circle (Example 7) 


Example 8. (Lead-lag Compensation) We give an example of lead-lag compensation design uti- 
lizing the developments described above. Let us consider the interval plant 

G(s) = 63s 3 + h? + bis + bo 

with coefficients bounded as follows: a 0 £ [5, 7], 63 E [.09, .11], 62 € [.9, 1.2], b\ E [. 8 , 1.5], bo E [.1, .3]. 
The objective of the design is to guarantee that the entire family of systems has a phase margin of 
at least 60° and a gain margin of at least 30dB. From Figure 28, we observe that the phase margin 
70° which is equal to the desired phase margin of 60° plus some safety factor can be obtained if the 
new gain-crossover frequency u>' c is at 0.35 [rad /sec]. This means that the phase-lag compensator 
must reduce the maximum magnitude of G(ju f c ) over the interval family to 0 [db]. Thus, we solve 

- 201 og 10 a = max \G(jw' c )\ = 28[db] 

G(ju/')€ G(ju' c ) 

and we have a = 0.0398. In order that the phase lag of the compensator does not affect the phase 
at the new gain-crossover frequency o>', we choose the value of 1 /aT to be at one decade below uj f c . 
Thus, 

T _ 10 _ 10 -7i7ft 7 q 

aw' c (0.0398)(0.35) 
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Therefore, the lag compensator obtained is 

1 + aTs _ 28.5714s + 1 
i(s) - 1 + Ta - 717 875a + 1 - 

We have now achieved approximately 70° of guaranteed phase margin and 25dB of guaranteed gain 
margin. To achieve the desired gain margin, we now wish to move the phase-crossover frequency lj" 
to 4.7 [rad/sec]. If the magnitude plot does not move, we can achieve the gain margin of 35db at 
this frequency. Thus, we solve 

— 10 log 10 a = -(35 - 25) = -10[db] 

and we have a = 10. Then, 

1 = = ^10(4.7) = 14.86 

and T = 0.0673. Therefore, the cascaded lead compensator is 

1 1 + aTs 1 (10)(0.0673)s + 1 _ s + 1.485 
“ a 1 +Ts “ 10 0.0673s + 1 _ s+ 14.859' 
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Prom Figure 28, we verify that the compensated system provides approximately 105° of guaranteed 
phase margin and 50dB of guaranteed gain margin. Therefore, the controller 


C(s) = C 2 (s)C l (s) = 


s + 1.5 28.5714s + 1 
s + 15 717.682s + 1 



-250 -200 -150 -100 -50 0 

Degree 

Figure 30: Nichols templates (Example 8) 


attains the design specifications robustly. Figures 29, 30 and 28 show the Nyquist, Nichols and Bode 
envelopes of the uncompensated and compensated systems. 

Example 9. (PI Compensation) Consider the plant 

G( s _ + l 

.02s 4 + d 3 s 3 + d-is 2 + .04s 
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with its coefficients bounded by given intervals as follows: 

Ti\ G [.35, .45], d 3 6 [.25, .35], b 2 G [.9, 1.1] 

The objective of the design is to guarantee that the entire family of closed loop systems has a phase 
margin of at least 45°. 

Here we design a PI compensator of the form 


C(s) = K P + — . 

s 

From Figure 32 we see that the new gain-crossover frequency u>' c should be moved to 0.034 [rad/sec]. 
Since the maximum magnitude at cj' c is 55 [db], we let 

G P (M) := max \G(ja>' e )\ = -201og 10 K v = 54.89[db] 


from which 


K p _ jQ-|G p (ju;^)![db]/20 _ ^q-54.89/20 _ Q.0018. 


In order that the phase lag of the PI compensator not effect the phase of the compensated system 
at a;', we choose the corner frequency Ki/Kp to be one decade below w* c - Thus, 


Ki = ^K P = ^^(0.0018) = 6.12 x 10 -6 . 


Therefore, the PI compensator obtained is 

, 6.12 x 10" 6 + .0018s 

c{a) = 7 • 

Figure 31 shows the Nyquist envelope while Figures 32 and 33 show the Bode and Nichols envelopes 
of the uncompensated and compensated systems. We can see that the phase margin specifications 
are robustly achieved. 



30 










4^8.2 Optimal Controller Parameter Selection 

In the design problem, we interchange the role of plant and controller, and consider a family of 
controllers stabilizing a fixed plant G(s). The design problem here is to select from the controller 
family the best parameter value according to some criterion. This type of problem is common in 
practice. Let us consider for example that this parameter selection is to be made to maximize gain 
margin or phase margin. The set of controller parameters may be given in terms of bounded values. 
In other words, we have a interval family of controllers stabilizing a fixed plant. To maximize the 
gain margin (or phase margin) over this set, we explore the boundary results given earlier. 

In Section 4.3 we showed that the minimum gain margin (or phase margin) occurs over the 
extremal set. In practice, the maximum gain margin (or phase margin) will also frequently occur 
over the same set. To see this, we argue as follows. The Nyquist boundary of the system cuts the 
real axis over the range A to B (see Figure 34 and Figure 35). 



Figure 34: A system that delivers the maximum gain margin 



Figure 35: A system that does not deliver the maximum gain margin 


The point A corresponds to the minimum gain margin of the family. The point B is a potential 
candidate for the maximum gain margin (gain margin = 1/OB). Corresponding to the point J3, there 
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exists a set Cb of admissible (parameters lie within the given intervals) controllers whose Nyquist 
plots pass through the point B. Suppose that one of these controllers is such that the Nyquist 
plot does not cut the real axis at any other point (see Figure 34). It is clear that this controller 
delivers the maximum gain margin (1 /OB) and is the optimally robust gain margin controller in 
the family. We can see that the maximum gain margin will be unrealizable only if each and every 
controller in Cb happens to cut the real axis at another point also (see Figure 35). Based on the 
above arguments, we suggest the following design procedure. First, check whether the given family 
of interval controllers stabilizes the plant u(s) by using the GKT. If so, we use this set of controllers. 
Otherwise, we determine the parametric stability margin around the nominal values of controller 
parameters and create the largest interval stabilizing controller family around this nominal. From 
this controller set we can find the controller parameters that provide the maximum gain or phase 
margin. This can be done by generating the Nyquist envelope of the family using the development 
discussed earlier. If the margin obtained is not satisfactory, we reset the controller parameters to the 
new nominal corresponding to the point B and create a new interval stabilizing family. We repeat 
this procedure until a satisfactory margin is achieved or the improvement of the margin is negligibly 
small. The set of stabilizing interval controllers can be determined by many different methods; for 
example the locus introduced by Tsypkin and Polyak may be used. Of course, there is no guarantee 
that a globally optimum design will be achieved by this method. However, a satisfactory robust 
design will often result. 

Example 1 0. Consider the feedback system with the plant transfer function 

G(s) . = 250 

1 ; * .025s 2 + 1.75s + 30 

and controller transfer function 

CW - ? 

The nominal values of controller parameters are a 0 = 5, fio = 2, 70 = 1. 

The phase margin of the given system with the nominal values of controller parameters is 
19.55°. The objective is to tune the three controller parameters so that the resulting system has a 
phase margin of approximately 45°. We first find the 1 2 parametric stability margin 77 around the 
nominal values of parameters (a 0 ,/3o,7o) by using the method proposed by Tsypkin and Polyak. 
From this value of 77 , which was obtained as ri = 1, we create the family of stabilizing controllers in 
the form of an interval transfer function Ci(s). Using this interval controller family, we create the 
corresponding Nyquist envelope which shows the maximum obtainable phase margin of 29.85°. Using 
the formulas developed previously, we select the parameter values (ai,/?i,7i), and consequently the 
controller C(s,a i,/?i>7i) € Ci(s1 that produces the maximum phase margin of 29.85°. Since the 
resulting controller C(s, a \ , (3\ , 71 ) does not satisfy the given requirement, we proceed with a second 
iteration. The parametric stability margin around (ai,/3i,7i) we found is again I. We now create 
a new interval family C 2 ($) of stabilizing controllers. Using the same procedure as before, we find 
the controller C 2 (s,a 2 ,/32, 72) € C 2 (s). This produces the maximum phase margin of 43.76°. In the 
table below we present the successive designs through several iterations and Figure 36 shows the 
Nyquist plot of the optimal system for each iteration. 

iteration parameter ranges parameter selected phase margin 
i {a, (aufiuli) 

1 ([4, 6], [1,3], [0,2]) (4] l] 2) 29.85° 

2 ([3, 5], [2, 4], [1,3]) (3,4,1) 43.76° 


4^8.3 Discrete Time Control Systems 

The frequency templates of discrete interval control systems can also be constructed similarly. We 
still use the basic geometric facts regarding the addition, multiplication and division of complex 
plane image sets discussed here; however, the main simplifying tool used in the case of continuous 
systems, the GKT, is no longer applicable here. We illustrate this by an example. 

Example 11. Consider the discrete-time feedback system shown in Figure 37. Let the interval plant 
G (z) and the controller C(z) be 


G(z) = 


z + a 

z 2 + / 3 z + 7 


and 


C(z) = 


— 0.2z - 0.35 
z 3 
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Figure 37: A discrete time feedback system (Example 11) 


where 

q <E [0.8, 1.2] := [a-,a+], 0 <E [1.8, 2.2] := [/?-,/?+], 7 6 [2.8, 3.2] := [ 7 ", 7 +]. 
Then we have the extremal subset Ge(^) consisting of the following 12 systems: 

^ . z + (Aa~ + (1 — A)a + ) _ ,, z + (\a~ + (1 - \)a+) 


r ( -\- z + ( Xa +(l~A)a + ) 
Ge ' w = — — 


G'e,(=) = 


z + (Aa~ + (1 - A)a + ) 


Ge & ( z ) = 
Ge 7 ( z ) = 


3 ' z 2 + 0~z + 7 + 

z + a~ 

2 2 + (A/3- + (l-A)/3+)2 + 7 - 

z + a + 

z 2 + (A/3 - + (1 - A )0 + )z + 7 ~ 


G E 8 (z) = ~2 


Geh {%) 2 


z 2 + 0~z + (A 7 ~ + (1 - A) 7 +) 

z + a + 

z 2 + 0~z + (A 7 _ + (1 - A) 7 +) 


G E2 (z) = - 
G e*(z) = ~ 

G E e ( z ) = ^2 

Ges i z ) ~ ~2 


Z 2 + 0 + Z + 7 “ 

z + (Aa~ + (1 — A)a + ) 
z 2 4- 0 + z + 7 + 

z + a~ 

z 2 + (A/J- + (1 - A )0+)z + 7 + 

z + a + 

z 2 + (\0~ + (1 - A )0 + )z + 7 + 


Ge 10 ( z ) = ~2 


G Ei 2 (z) = 


z 2 + (3 + z + (A 7 _ + (1 — A) 7 +) 
z + a+ 

z 2 + /3+z + (A 7 “ + (1 - A)t+) 


where in each case A ranges over [0, 1]. 

By searching over the family Ge(^) with z = e^ T for u) E [0, oo) and a fixed T, we can obtain 
the frequency templates of the system. For illustration, the Bode envelopes of the above discrete 
time feedback system are given in Figure 38 for T =1. The Nyquist and Nichols envelopes can also 
be generated similarly. 
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5 MULTILINEAR INTERVAL SYSTEMS 

Consider the feedback configuration shown in Figure 39. F(s) is a fixed controller and G\ ( s ) and 



Figure 39: Interconnected Feedback System 


G' 2 (s) are independent subsystems containing parameter uncertainty. Let 


F(s) := 


Fi(s) 

F2(S)’ 


G i(s) := 


Pni£t 

P2i(s)' 


G 2 (s) := 


P M 
F22W 


(36) 


If the subsystems &i(s) and G 2 (s) contain independent parameters it is reasonable to model them 
as interval or linear interval systems G;(s), i = 1,2. The characteristic polynomial of this feedback 
system is* 

' " 6(s) := F 1 (a)P n (a)PM + F 2 (s)P 21 (s)P 22 (s). (37) 

The open loop transfer function 


T°(s) = F(s)Gi(s)G 2 (s) — 


F\^]Pn{s)P^ 

F 2 (3)P 2 i(s)P22(s) 


and the closed loop transfer function 

F(s)Gi(s)G 2 (8) = F(s)P u (a)PM 

W l + F(s)G l (s)G 2 (s) F 1 (s)P 11 (s)P 12 (a) + F 2 (s)P 21 (s)P 22 (s)- 

These transfer functions have numerator and denominator polynomials that are multilinear functions 
of interval polynomials. In addition, the numerator and denominator in T°(s) have independent 
interval polynomials (parameters), but in T c (s) the numerator and denominator contain common 
interval polynomials. 

For this class of systems we are interested in the following types of questions: 

1) Does F(s) robustly stabilize the system or not? 

2) If F(s) does robustly stabilize the system, what are the worst case gain margin, phase margin, 

parametric stability margin, stability margin and performance measured in terms of JToo 

norms, as the parameters range over the uncertainty set? 

3) How can one construct the Bode magnitude and phase, and Nyquist plots of various transfer 

function sets such as T°(s) and T c (s ) generated by the uncertain parameters? 

These questions were addressed in the previous sedction for an arbitrary stability region using the 
Mapping Theorem as a computational tool. In this section, we focus on the case of Hurwitz stability 
and show how the Generalized Kharitonov Theorem (GKT) can be extended to this multilineax case 
to provide a great deal of simplification and computational efficiency. We first derive a multilinear 
version of GKT. This generalization provides us with an extremal test set for multilinear interval sys- 
tems, with drastically reduced dimension of the parameter space, from which all the above questions 
can be answered. It will tell us, essentially that as long as the dependencies are multilinear, worst 
case stability margins and performance can be calculated for an arbitrary control system containing 
interval subsystems by replacing each interval subsystem G*(s) by the corresponding extremal set 
of systems G^s). The polytopic approximation derived in the last section can then be used on this 
extremal set to give a highly efficient computational solution to the design and analysis questions 
posed above. In the next section we describe the extension of GKT to the multilinear case. 
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5.1 MULTILINEAR INTERVAL POLYNOMIALS 

To avoid notational complexity, we first consider the simplest multilinear polynomial form motivated 
by (37). We consider the Hurwitz stability of the characteristic polynomial family of the form 

S(s) := FiWPuWPiaW + F 2 (s)P 2l (s)P 22 (s) (38) 

where Fi(s) are fixed real polynomials and Pij(s) are real interval polynomials with independently 
varying parameters. Let p denote the ordered set of coefficients of the polynomials {Pn(s) y P\ 2 {s)P 2 \ , (s)P 22 (s)}. 
We assume that each coefficient varies in an independent interval, or equivalently that p varies in 
an axis-parallel box II. The dimension of the parameter space here is equal to the number of 
independently varying parameters contained in these polynomials. We let Pij(s) denote the in- 
terval polynomial family and ICij(s) and Sij(s) denote the respective Kharitonov polynomials and 
Kharitonov segments. The family of uncertain polynomials is represented as 

A(s) := F\ (s)P ii (s)P \ 2 (s) + Fi(s)P 21 (s)P 22 (s). (39) 


Define 

Ae(s) := Fi(3)<S u (s)5 12 (s)+F 2 ( S )K: 21 (s)X: 2 2(s), A|(s) := F 1 (s)ICn(s)ICMs)+F 2 (s)S 21 (s)S22(s) 

and introduce the extremal manifolds 

Ae(s) := Ag(s) U A|(3). (40) 

Lemma 2. Under the assumption that every polynomial in Afs) is of the same degree and the 
parameters (coefficients) in the interval polynomials Pn(s), Pi 2 (s), P 2 i(s), P 22 (s) are independent , 
Ms) is Hurwitz stable if and only if Ae(s) is Hurwitz stable . 

For the general case, we have the multilinear interval polynomial 


M*) = F x (s)Pn(s) • • • Pi ri («) + ■■• + F m (s) P mX (a) • • ■ P mrm (a) (41) 

where Fi(s) are fixed and Pij{s) are interval polynomial families. The arguments used in the earlier 
case carry over to this case and the result will be presented without repeating the proof. We define 

Ag (3) := F\ (s)/Cn (3) * • • fci ri (3) H h P/-i (s)fCi-\ t i (3) • • ■ (s) 

+F/(s)«S/,i (3) • • • S/, r , (3) + F/+1/C/+14 (s) * • * A^f+i,r/ + i + 1 - F m (s)fC mt 1 (s) * • • £ m ,r m (5) ( 42 ) 

and introduce the extremal manifolds 


A E ( S ):=U, m =1 A'(»). (43) 

Theorem 11. Under the assumption that every polynomial in A ( 3 ) is of the same degree and the 
parameters (coefficients) in the interval polynomials Pij(s) are independent, A(s) is Hurwitz stable 
if and only if Ae (3) is Hurwitz stable. 

Remark 4- The theorem is also valid when some of the polynomials F*(s) are complex or quasipoly- 
nomials This follows from the corresponding fact that GKT holds when some of the fixed polynomials 
Fi(s) are replaced by complex polynomials or quasipolynomials. 

The next example illustrates the construction of the set of manifolds. 

Example 12. (Extremal Manifolds for Multilinear Systems) Consider the characteristic poly- 
nomial 

S(s) = F 1 (s)P n (s)Pi 2 (s) + F 2 (s)P 2 i(3)P 22 (3) 
where the fixed polynomials Pi (3) are 

Fi(s) = s + / 1 , F 2 (s) = s + f 2 


and the interval polynomials are 

P11 ( 5 ) — a 2S 2 +ois + oq, Pi 2 (s) = b 2 s 2 + 6 iS + &o> P 2 i(s) = c 2 3 2 +cis+co, P 22 (s) = d 2 s 2 + d\s + do, 
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with all coefficients varying independently. The Kharitonov polynomials corresponding to I\ i (a) , 
Pi2(s), P2i{s) and P 2 2(«) are: 


Kp n (s) = a^s 2 +a x s + o 0 , 
K p n (s ) = a^s 2 + a+s + aj, 
k p 12 («) =b^s 2 + fcj - s + fcj, 
K P 2l ( s ) = c 2 s2 + C t S + C 0 ' 

Kp 22 (s) = d£s 2 + d 1 s + d 0 , 
Kp 22 ( s ) = d 2 s 2 + d^s + . 


K 2 Pll (s) = a^s 2 + a^s + a 0 , 
Kp 12 (s) = b%s 2 + 6fs + 6o , 
Kp l2 ( s ) = b 2 s 2 + s + b$, 
k p 21 (*) = Cj s 2 + cf a + a£, 

Kp 2 3 (a) = ^ s 2 + df 3 + d 0 , 


Kp u ( s ) = 02 s 2 + of s + aj , 
^p 12 ( s ) = f>ts 2 + bfs + b^ , 
k p 2x ( s ) = 4 s 2 + cf s + c<7, 
K ki («) = c 2 s2 + C 1 S + C 0 - 

■K p 23 (5) = d*2 & T dj s ■+• dg » 


The sets of segments joining appropriate pairs of Kharitonov polynomials can also be obtained. Here 
we give only the segments corresponding to the interval polynomial Pi 1 (s); others can be similarly 
obtained. 

Sp u (s) = *K l pJs) + (1 - A )K 2 Pn (s), S 2 Pll (s) = XK'pJs) + (1 - A )K%Js), 

S Pll (s) = A K 2 Pll (s) + (1 - A)-Kp n (s), S 4 Pu (s) = A K 3 Pll (s) + (1 - X)K 4 P Js). 

From the above, we now can write the manifolds: 


Ax (a) = {F 1 (s)K i Pll (a)K } P Js) + F 2 ( S )S k p 3 l (s)S l p 3 ,(s) : (i,j,k,l) €4x4x4x4} 
A 2 (s) = {F^S^ (s)S 2 p Js) + F 2 (s)Kp 21 ( s)K l P ^s ) : (i, j, k,l) € 4 x 4 x 4 x 4} 


and 

Ae(s) = Ai U A 2 . 

As we can see, the total parameter space is of dimension 12. However, our problem is now reduced 
to checking the stability of 512 two-dimensional manifolds. Notice that each manifold remains of 
dimension two, even though the dimension of the parameter space can be increased arbitrarily. The 
dimension of the manifolds is increased only if the number of interval polynomials in a product term 
is increased, 

5.1.1 Dependencies Between the Perturbations 

The theorem is stated assuming that the polynomials Pa (s) perturb independently. In an intercon- 
nected multiloop control system, it will in general happen that some of the polynomials Pij(s) are in 
fact identical (see, for example, Exercise 11.1). Such dependencies can be easily handled. To avoid 
introducing cumbersome notation we illustrate the procedure with an example. 

Consider the following multilinear interval polynomial family with dependencies between the 
perturbations: 

6(s, p) := F 1 (s)Pn(s)P 12 (s) + F 2 (s)P 21 (s)P 22 (s) + F 3 (s)P 3l (s)P 32 (s) 

where P u (s) = P21W and P 22 (s) = ^32(5) and each polynomial Pij(s) is interval. If we rewrite 
embedding the above constraints, we have 

6 (s,p) = (44) 

Let us first fix Pn(s) and P 3 i(s) and apply GKT. This tells us that the Hurwitz stability of the set 
A (s) is equivalent to stability of the sets 

Ii (s) = {F 1 (s)Pn(s)S l2 (s) + (F 2 (s)P u (s) + F 3 (s)P 3l (s)))C 22 (s)} 
h(s) = {^WPnW/CiaW + (F 2 (8)Pn(8)+F 3 (s)P 31 (8))S 22 (8)} 

for each (Pn(s),P 3 i(s)) E Pu(s) x P 3 i(s). We now apply GKT again to each of the sets Ii (s) and 
I 2 (s) letting Pn(s) and P 3 x(s) now vary. This leads to the condition that the robust stability of 
A(s) is equivalent to robust stability of the manifolds: 

Ai(s) = {6(s,p) : Pn E 4S11, P12 € $12, P22 € ^22, P31 € /C 3 i} 

^2(5) = {£(s,p) : Pn E Siu Pi 2 € )Ci2i P22 € <$22* P31 ^ ^31} 

A 3 (s) = {J(s,p) : Pn E 1 C n, P12 € <Si2> P22 € ^221 P31 € 5 3 i} 

A4 (s) = {<5(s,p) : Pn € fCiu P12 € fC\2, P22 ^ <$22, ^ 31 ^ ^ 31 ) • 
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Since Aa(s) is a polytope, we now can apply the Edge Theorem to conclude that this is stable if 
and only if its exposed edges are. These exposed edges are 

{£(s,p) : P \ i G /Cn, P\2 G S 12 , T 22 € K, 22, P31 £ £ 31 } (45) 

{<5(s,p) : Pxx G £n, Pi 2 G £12, T22 6 /C22, ^31 € £31 } . (46) 

It is easy to see that the manifolds (45) and (46) are contained in Ai(s) and A 4 (s) respectively. 
Therefore the set of manifolds that finally need to be checked is 

AbW := A x (s) U A 2 (a) U A 4 (s). (47) 


5.1.2 Interval Polynomial Matrix 

An important special case of dependent perturbations is that of interval polynomial matrices. Let 
M{s) be an n x n polynomial matrix whose ij th entry is the polynomial Pij (s). The characteristic 
polynomial associated with this matrix is det[M(s)]. We assume that each Pij(s) belongs to an 
interval family P^(s) and the parameters of each (s) are independent of all others. Let M(s) 
denote the corresponding set of matrices. We will say that the set is stable if det[M(s)] is Hurwitz 
stable for every M(s) e M(s). To test Hurwitz stability of the family of characteristic polynomials 
so obtained we have the following result. Let T be the set of n x n permutation matrices obtained 
from the identity matrix. Corresponding to each T G T, introduce the set of polynomial matrices 
Mt(s) where Pij{s) ranges over JCij(s) if the ij th entry of T is 0 and Pij(s) ranges over Sij(s) if the 
ij ih entry of T is 1. Let M*(s) denote the collection of matrices M-r(s) obtained by letting T range 
over all permutation matrices T. 

Theorem 12. M(s) is Hurwitz stable if and only i/M*($) is Hurwitz stable. 

5.1.3 Boundary Generating Property of the Extremal Manifolds 

We have established that Hurwitz stability of the family A(s) is equivalent to the stability of the 
manifolds Ae(s). This equivalence also follows from the following boundary result relating the 
image sets A (ju;) and AeC?^). 

Theorem 13. 

dA(ju) C A e (jv) (48) 

The boundary result established above implies that the problem of checking robust stability is now 
reduced to verifying that the origin is excluded from the set AeO’uO for each u in [0, 00 ). This is 
a great deal simpler than the original problem because of the greatly reduced dimensionality of the 
set AeW relative to the original set A (s). However, this verification is still not easy because the 
set Ae(4 is multilinear in the parameters A*. At this point, the Mapping Theorem can be brought 
in and used to approximate Ae (jw). Indeed, since Ae(s) depends multilinearly on the parameters 
A ij associated with the Kharitonov segments, the vertex set of A (s) can be generated by setting the 
Pij(s) to the corresponding Kharitonov polynomials. Let 

A k (s) = {<5(«,p) : Pij(s) = Kly(s), i 0 €4, i € m, j € »Vn}- ( 49 ) 

We can also introduce the poly topic set consisting of convex combinations of the vertex polynomials: 
A(s) = {At><(«) + (1 - A)Uj(s) : Vi(s),Vj(s) e A K (»), A €[0,1]}. (50) 

It follows from the Mapping Theorem now, that 

co AeO'u;) = co A K (jo;) = A K (i^). 

Therefore, the condition that A E (j^) excludes the origin can be replaced by the sufficient condition 
that Ak(.?^) exclude the origin. Since this latter set is a polytope, this condition can be verified by 
checking that the angle subtended by the set AkO’cj) at the origin, is less than * radians. 

Therefore we have proved the following result. 

Theorem 14. The family A (s) is Hurwitz stable if it contains at least one stable polynomial and 
satisfies 

1)0$. co A (ju) for some lj 2) $a k (J u ) < f° r a h w G [0, oo). 
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This result states that the Hurwitz stability of A (5) can be determined by checking the phase 
difference of the vertex polynomials corresponding to the Kharitonov polynomials along the ju) axis. 
This is useful in view of the fact that the number of Kharitonov vertices is fixed whereas the vertices 
of II increase exponentially with the dimension of the parameter space. We illustrate these results 
with examples. 

Example 13. Consider the interconnected feedback system shown in Figure 40. 



Figure 40: Interconnected Feedback System (Example 13) 


Let 


F(s) 

Q(s) 


Fi($) _ s + 2 jMg) _ s 2 + s + 1 

F 2 (s) s + 1’ 8 P 2 (s) s 3 + a 2 s 2 + 4s + a 0 ’ 

Qi(s) _ 6.6s 3 + 13.5s 2 + 15.5s + 20.4 
Q 2 (s) s 3 + bzs 2 + 3.5s + 2.4 


and let the set of parameters p = [a 2 ,a 0 ,& 2 ] vary as follows: 
a 2 G fa ,aj] = [-3.625,-2.375], a 0 € fa,a+] = [1.375,2.625], b 2 G fa fa] = [2.875,4.125]. 


The characteristic polynomial of the system is: 


<5(s,p) = s 7 + (7.6 + q, 2 4* 5 2 )s 6 + (40.8 + Q 2 + b 2 -\- u 2 & 2 )s 5 -I- (85.7 + oo + 3.5o 2 + 46 2 -+■ a 2 5 2 )s 4 
+(137 + do + 5.9a 2 + 46 2 + cio5 2 )5 3 + (158.3 + 3.5ao + 2.4d 2 + 00 ^ 2 )^^ 

+ (101.8 + 5.9a 0 )s + (40.8 + 2.4a 0 ). 


We verify that the following polynomial in the family is Hurwitz: 


S(s, p = [-3, 2, 3.5]) = s 7 + 8.1s 6 + 30.8s 5 + 80.7s 4 + 142.3s 3 + 165.1s 2 + 113.6s + 45.6 


The parameter sets corresponding to the Kharitonov polynomials are 


{ ( a 2 ) a 0 5 ^2") > ( a 2 > a J 1^2 ) » ( a 2 » a o > ^2 ) > ( a 2 i a 0 ’M} • 

The set of Kharitonov vertex polynomials is 

A k («) = {^/Ci(5),^ 2 (s),rfic 3 (s),<5K4(5)} 

where 

6k 1 (s) = s 7 + 9.35s 6 + 32.7531s 5 + 85.4656s 4 + 146.5344s 3 + 163.0844s 2 + 109.9125s + 44.1 

6 K2 (s) = s 7 + 8.10s 6 + 34.4719s 5 + 83.4344s 4 + 139.8156s 3 + 161.3656s 2 + 109.9125s + 44.1 

<5*3 (s) = s 7 + 8.10s 6 + 26.3469s 5 + 77.1844s 4 + 145.5656s 3 + 169.6156s 2 + 117.2875s + 47.1 

<5* 4 (s) = s 7 + 6.85s 6 + 29.6281s 5 + 76.7156s 4 + 137.2844s 3 + 166.3344s 2 + 117.2875s + 47.1 

One approach is to check the Hurwitz stability of all convex combinations of these polynomials. This 
in turn can be done by using the Segment Lemma. Alternatively, we may check the phase differences 
of these vertex polynomials. If the maximum phase difference is less than 180° (ir radians) for all 
w, the origin is excluded from convex hull of Ag(j^) for all u. We show the convex hulls of image 
sets in Figure 41 for illustration. 

We confirm robust stability by verifying (see Figure 42) the maximum phase difference never 
reaches 180° for all u. We also note that if we had applied the Mapping Theorem directly to the 
three dimensional parameter space p, we would have had to check 8 vertices as opposed to the 4 
vertices checked here. This reduction is due to the application of the multilinear version of the 
Generalized Kharitonov Theorem. 
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Figure 41: Convex hulls of image sets (Example 13) 



Figure 42: $a k (w) vs. ^ (Example 13) 


5.2 PARAMETRIC STABILITY MARGIN 

Consider again the family of polynomials 


A (s) = F x (s)Pn (s) P lri (*) + ••■ + Fm(s) P m i (s) ■ • • P mrm (s) . (51) 

Let p ij denote the vector of coefficients of the polynomials W- Each such coefficient belongs 
to a given interval and the corresponding box of parameters is denoted by 11^. Then the box of 
uncertain parameters is 

rr := iin x rii 2 x ••• x n mrm . (52) 

The extremal manifolds Ae(s) are defined as follows: 


Ae( 5 ) := F\(s)/Cn(s) * * * 1C i n ($) + ■*• + ( s ) 

+Fi(s)Si y i(s) • • *<Sf, n (s) + ( 5 ) • * * /Cj+i >n+1 + • * • 4- E m (s)/C m> i ( 5 ) * * • £ m , rm (s) (53) 

A e (s) := Ag(s). (54) 
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( 55 ) 


The parameter space subsets corresponding to A' e (s) and A E (s) are denoted by Ilf and 

m 

n E := (Jn,. 

i=l 

Let IIk denote the parameter vector set corresponding to the case where each polynomial Pa (s) set 
to a corresponding Kharitonov polynomial. We write 

A(s) = {<5(s, p) : pen}, A E (a) = {<*(«, p) : p € n E }, A K (s) = {<5(s, p) : p € n K }. 

In this section, we show that for the family A (s) in (51) the worst case parametric stability 
margin over the uncertainty box II occurs, in fact, on the set IIe. 

Let || • • * || denote any norm in R n and let V u denote the set of points u in R n for which 6(s> u) 
is unstable or loses degree (relative to its degree over II). Let 

p{ p) = inf Up - u||p 

denote the radius of the stability ball (measured in the norm || • ||) and centered at the point p. This 
number serves as the stability margin associated with the point p. If the box II is stable we can 
associate a stability margin with each point in II. A natural question to ask is: What is the worst 
case stability margin in the norm || • || as p ranges over II? The answer to that question is provided 
in the following theorem. 

Theorem 15. (Extremal Parametric Stability Margin) 

inf/>(p)= inf p(p)- (56) 

pen penE 

Example 14 ■ Consider the system in Example 13 with nominal values 

a% = -3, a° 0 = 2, b° 2 = 3.5. 

We compute the maximum parametric stability margin e* around the nominal values as follows: 

a 2 € [ fl 2 ~~ £, + c], do E [oq ~~ a o T e], 62 ^ [ 5 ^ — e, &2 T- e]* 

Prom the four vertices given in Example 13, we have six segments bounding the convex hull of the 
images. To check the stability of these segments, we apply the Segment Lemma with incremental 
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Figure 43: Image set for e = 0.63 (Example 14) 

steps of c. This gives e * = 0.63. The image set of the characteristic polynomial with e* is shown in 
Figure 43. The figure shows that the image set is almost touching the origin and thus e* = 0.63 is 
the parametric stability margin. 
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5.3 MULTILINEAR INTERVAL SYSTEMS 

In the rest of this section, we will be dealing with transfer functions containing interval parameters 
with multilinear dependency. The transfer function in question could be embedded in a feedback 
control system and we will be interested in determining robust stability as well as worst case stability 
margins and performance measures of such systems. We begin by first focussing on a multilinear 
interval system , namely one whose transfer function is a ratio of multilinear polynomials with in- 
dependent parameters. To be specific we will consider single-input, single-output, proper, stable 
systems with transfer function of the form 


G(s) 


7(f) 

Sis)' 


Here 

7(5) = Hi(s)Ln (s)Li2(s) • • * L\ ri (s) + + H Tn (s)L m i(s)L m 2{s ) • • • L mrrn ( s ) 

where the polynomials Hi(s) are fixed and the polynomials Lij(s) E Lij($) are independent real 
interval polynomials. Let 1 denote the ordered set of coefficients of the interval polynomials; 1 varies 
in a prescribed axis parallel box A; the corresponding family of polynomials 7(3) is denoted by r(s). 
Similarly we suppose that 

S(s) = Fi (s)P n (s)Pi2 (s) • * • Pi ri (s) + + F m (s)P ml (s)P m2 (s) • • • P mfm (s) 

where the polynomials Fi(s ) are fixed, the polynomials Pij(s) are real interval polynomials with 
the vector of coefficients denoted by p varying in the prescribed box II. The resulting family of 
polynomials <5(s) is denoted A (5). We also denote explicitly the dependence of <5(s) on p and of v(s) 
on 1 by writing <5(s, p) and 1/(5, 1) whenever necessary. We make the standing assumption. 

Assumption 1. 

Al) Parameters p and 1 are independent. 

A2) 7 ( 5 , 1 ) and (5(5, p) are coprime over (p, 1) E II x A. 

A3) 6(jcj, p) # 0 for all p E II and each u > 0. 

Later in this section we show how to deal with the situations when assumption Al does not hold. 
To display the dependence of a typical element G(s) of G on 1 and p we write it as G{s, p, 1): 


G(5,p,l) 


tOM) 

<5(s,p) 


(57) 


We form the parametrized family of transfer functions 

G(«) = {<?(«, P.1) : (p,l)e(nxA)} = ^. 


(58) 


In order to apply frequency domain methods of analysis and design to the family of systems 
G(s) it is necessary to obtain the image set G(ju>). We first show how the boundary of this set can 
be evaluated. We proceed as in the linear case by determining an extremal multilinear interval family 
of systems Ge(s) . Introduce the Kharitonov polynomials and segments associated with the Pii(s) 
and Lij(s) respectively, and construct the extremal polynomial manifolds Ae(«) and Te(s) and the 
vertex sets A K («) and r K (s) as in Section 5.1. Let IIe, Ae and IIk, Ar denote the corresponding 
manifolds and vertices in II and A respectively: 

r E («) = : 1 G A e } , T K (s) = ( 7 («,1) : 1 e A K ) 

A e (s) = {<5(s,p) : P € n E } , A k (s) = {6(«,p) : P € n K } . 

The extremal set Ge(s) is then defined as 

tM . 

<S(«,P) ’ 

Using our compact notational convention we can write 



(1 E Ak, p € IIe) or (I E Ae, p € IIk) j ■ 


(59) 


G e (s) = 
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Theorem 16. Under the Assumption 1, 


dG(juj) C G E (jv) 


for all lj € [0, oo). 

Now suppose that G(s) is part of the control system shown in Figure 44. 



Figure 44: A unity feedback system 


Define transfer functions: 



U{S) - F(s) 

e(s)~ HSh 


(60) 

T°(s) := ^ = F(s)G(s), 

e[s) 

T e (s) - = 

T (S) - r(s) 

1 

(61) 

1 + F(s)G(s) ’ 

_ u(s) F(s) 

t«( 3 ) — — 

F(s)G(s) 

(62) 

’’ r(s) l + F(s)G(s)’ 

T (3) - r(s) ~ 

1 + F(s)G(s) ' 


As G(s) ranges over the uncertainty set G(s), these transfer functions vary in corresponding sets. 


T°(s) := {F(s)G(s) : G(s) G G(s)}, 

F(s) 


n.):={ T 


+ F(s)G(s) 


G{s) G G{s) 


T '<»> : = {ttf5 

T» :=| T 


F(s)G(s) 
F(s)G{s) 

+ F(s)G(s ) 


: G(*)eG(*)j 

: G(s ) G G(s)l (63) 


It turns out that the boundary of the image set, at s = ju, of each of the above sets is generated by 
the extremal set Ge(s). Introduce the extremal subsets: 

T°e(°) ■■= {F(s)G(s) : G(s) € G E (s)}, T| (s) := j j— : G(s) G G E (a)J(p4) 

TSW = { i ■ ° w e Ge(,) } ■ T|w := IS : GW e GE(>) } (li5) 


Theorem IT. For every to > 0, 

(a) dT°(ju) C T° E (ju) (b)dT*(ju>)cT%(ju>) (c) dT*(ju>) C Tgtfw) (d) BWiju) C T v E (ju>) 

All the boundary results related to the Nyquist and Bode envelopes given in the previous section 
carry over to the multilinear case with the corresponding extremal set Ge(s). In the following 
subsection we show that these boundary results hold for a much larger class of functions of G(s). 


5.3.1 Extensions of Boundary Results 

In defining our multilinear interval system we had assumed that the numerator and denominator 
parameters of G(s) are independent. However for closed loop transfer functions this assumption of 
independence does not hold. Nevertheless the boundary generating property of the set Ge(s) still 
holds for these closed loop transfer functions. It is natural to attempt to generalize this bound- 
ary generating property to a large class of functions where dependencies between numerator and 
denominator parameters occur. Such dependencies invariably occur in the transfer functions associ- 
ated with multiloop control systems. We begin with a multilinear function Q(s) of several transfer 
functions £*($), i = 1,2 •••,<?. Let us assume that each G‘(s) itself lies in a multilinear interval 
family of systems G *(s) defined as in this section, with independent parameters in the numerator 
and denominator. We also assume that the parameters in G l (s) are independent of those in G J (s), 
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i ^ j. Note that if we regard Q(s) as a rational function its numerator and denominator polynomials 
contain common interval parameters. Let 

Q (s) := {Q(s) : &(s) £ G*(s), i = 1,2 • <7}. 

We wish to determine the complex plane image set of the family Q(s) evaluated at s = jco : 

Q (ju>) := {Q(jv) : (?<(*) E G i (s), i = 1,2 ■■•*}. 

Let Gk(s) denote the extremal subset of G‘(s) and introduce 

Qe(M := {Q(j^) : G*(s) € G l E (s), i = 1,2 * 

Then we can state the following boundary result. 


Theorem 18. 

dQ (ju) C Qe(jw)- ( 66) 

We can generalize the above property even further. Let Q(s ) be as above and consider a linear 
fractional transformation (LFT) T(Q(s)) defined by arbitrary functions j4(s), B(s),C(s) and D(s) 


T(Q(s)):= Qi(s) := 


A(s)Q(s) + B(s) 
C(s)Q(s) + D(s)' 


Introduce the set 


T(Q(s)) := { 


yl(s)g(s) + B(s) 
C{s)Q(s) + D(s) 


Let us now impose the restriction 


Q(s) e Q(s)| . 


A(ju)D(ju) - B(ju)C(ju) / 0. 


(67) 


Under the above restriction T(Q(ju)) is a LFT of Q(juj) and thus carries boundaries onto boundaries. 
Thus we know that the boundary of Qi(jcj) is also generated by the extremal systems G^(s). 


Theorem 19. Let Q(s) be as defined above and let T(Q(s)) be an LFT. Then, assuming (67) holds 

dT(Q(ju)) C T(QeCM). 


Remark 5. We remark that any LFT of £?i(s) as well as sums and products of LFT’s continue to 
enjoy the boundary generating property of the extremal systems. In this way a large class of transfer 
functions occurring in closed loop systems can be handled. 


5.3.2 Computation of G E (ju>) 

As we have seen, the determination of the set G (ju;) reduces to evaluating Ge(j^), a set of smaller 
dimension. Nevertheless, Ge{Ju>) is still a multilinear function of the uncertain segment parameters. 
The only way to exactly evaluate this set is by gridding over the uncertainty box. In general this 
procedure is computationally expensive or even infeasible. Fortunately, as we have seen earlier, 
the “concave” property of the image set of multilinear interval polynomials given by the Mapping 
Theorem allows us to overbound G E (j^) by a ratio of unions of convex polygons. Recall that 

<68> 

From the Mapping Theorem we have 

T e (jw) c CO A E (ju;) C CO Ak (ju). 


Now consider G E (j'u;) defined by replacing r E (jt<;) and A E (jci;) by co T kU^) and co Ak ( jcj) 
respectively. In other words, 


Ge(jcj) 


-( 


w \ 

co Ak (ju)J 


U 


' co r K (ju>) \ 

, AkO'w) ) 


(69) 
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It is clear that Ge(j^) overbounds Ge(jo;): 

G e (j M C G e Uu>). (70) 

The evaluation of Ge(j^) is relatively easy because it consists of a union of one-parameter families 
of transfer functions of the types 

XUiiju) + (l - \)U 2 (ju>) U{ju>) 

V(ju>) ° r AVi CM + (1 - A)V 2 (jlj) ■ 

The union of these one-parameter families gives rise to Ge ( ju;) which overbounds the boundary of 
G(jcj). The tightness of the approximation can be improved as we have seen before, by introducing 
additional vertices in the parameter set II x A. The guaranteed gain and phase margins of a control 
system containing parameter uncertainty may be obtained from these overbounded sets as shown in 
Figure 45. The image at s = jcv of each of the transfer function sets associated with the feedback 






(c) (d) 

Figure 45; Guaranteed gain and phase margin using overbounded sets 

system considered in Figure 44 can be overbounded by replacing G(juj) by Ge in fact we 
can do the same for any linear fractional transformation of G(s). An example of this calculation 
follows. 
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Example 15. (Nyquist, Bode and Nichols Envelopes of Multilinear Systems) The purpose 
of this example to show how to construct the frequency domain envelopes (Nyquist, Bode and Nichols 
envelopes) of the multilinear interval family. Let us recall the system given in Example 12: 


(s 4- 2)(s 2 + £ + l)(6.6s 3 + 13.5s 2 + 15.5s 4 20.4) 
n — — l' (s 4- l)(s 3 + a 2 s 2 +4 s 4- a 0 )(s 3 + & 2 s 2 -f 3.5s -f 2.4) 

G(s) 


where 

a 2 € [a^ ,a+] = [-3.625,-2.375], a 0 <E [o 0 ",a 0 + ] = [1.375,2.625], & 2 € [hj.fcj] = [2.875,4.125]. 


From Theorem 16, it is enough to consider GeO’u;). Since this particular example has parameters 
only in the denominator, we have 


d(F(jw)G(j«>)) C F(jw)G b (jw) = 


(s + 2)(s 2 4- s + l)(6.6s 3 4- 13.5s 2 4- 15.5s + 20.4)|$^ 
AeO'o;) 


where Ae (ju) consists of the six line segments joining the following four vertices: 

(s 4 l)(s 3 4 a^s 2 4 4s 4 a,Q ) (s 3 4 b^s 2 4 3.5s -1- 2.4)1$— 

(s -f l)(s 3 4- aj s 2 4 4s 4- )(s 3 4- b^s 2 4- 3.5s 4- 2.4)|$ = j u; , 

(s 4" l)(s 3 4 (I 2 s 2 4 4^ 4 Uq” )(s 3 4 b^ s 2 4- 3.5s 4- 2.4) \ 8= j u)f 

(s 4- l)(s 3 4- s 2 4 4s 4- aj^)(s 3 4 5^s 2 4 3.5s 4- 2.4)|*=j W . 

Figures 46, 47, and 48 are obtained accordingly. 



Figure 46: Nyquist envelope (Example 15) 


5.4 Hoo STABILITY MARGIN 

In this section we use the boundary results derived above to efficiently deal with many frequency 
domain measures of performance and, in particular, to determine worst case performance over the 
parameter set associated with a multilinear interval systems. We will use the standard notation: 
C+ := {s € C : Re(s) > 0}, and i?oo(C+) will represent the space of functions f(s) that are 
bounded and analytic in C+ with the standard H norm, 

ll/lloo = sup |/0’w)|. 
u6R 
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Let us consider the multilinear interval family of systems G(s) defined earlier (see (58)) and let 
us assume that the entire family is stable. To determine the unstructured stability margin of 
control systems containing the family G(s) we need to determine the supremum of the H 0 0 norm of 
certain transfer functions over G(s). Since the Hqo norm-bounded perturbations provide uniform 
perturbations at all frequencies, it is desirable to shape this perturbation by introducing a weight. 
Let the weight W(s) be a scalar stable proper transfer function 


W(s) = 


n w (s) 

(^) 


To start with, let us look at two specific robust stability problems involving mixed parametric- 
unstructured uncertainty: 

Problem I: Consider the configuration in Figure 49, where W(s) is a stable proper weight, G(s) is 

a stable multilinear interval family of systems, and A P is any H ^ perturbation that satisfies 
||AP|| < a. Find necessary and sufficient conditions for stability of the family of closed loop 
systems. 



Figure 49: 


Problem II: Consider the feedback configuration in Figure 50, where VF(s) is a stable proper 

weight, A P is any perturbation that satisfies I |AP|[ < a, and C(s) is a controller that si- 

multaneously stabilizes every element in the set G(s). Find necessary and sufficient conditions 
for stability of the family of closed loop systems. 



Figure 50: 


The above problems axe generalized versions of standard Poo robust stability problems where 
a fixed plant is considered. Here, the worst case solution is sought over the parameter set II x A. 
We solve this problem by using the Small Gain Theorem and determining the worst case over the 
uncertainty set. The solution is accomplished by showing that the J?oo norms in question attain 
their supremum value over the extremal set of transfer functions Ge(&) C G(s) defined in (59). We 
can now state the main result of this section. 

Theorem 2(h (Unstructured Stability Margins) 

1) The configuration of (Figure 49) will be stable if and only if a satisfies 

1 _ . 

a “ su P»6G e Halloo a °’ 
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2) The configuration of Problem II (Figure 50) will be stable if and only if a satisfies 

1 * 

Q - sup 96Ge 11^(1+ pC)-loo :-<V 

Remark 6. The quantities a* and a* serve as unstructured Hqo stability margins for the respective 
open and closed loop parametrized systems treated in Problems I and II. 

Remark 7. In practice we can further replace Ge(s) by the image set overbounding poly topic family 
Ge(s), and obtain lower bounds on these margins from Ge(s). We can also obtain upper bounds on 
these margins from the set of Kharitonov vertex systems of the family. Furthermore, the extremal 
results stated above also hold for any LFT of G(s). 

5.5 NONLINEAR SECTOR BOUNDED STABILITY MARGIN 

We now consider the effect of nonlinear perturbations on the multilinear interval family G (s) defined 
in the last section (see (58)). The nonlinear perturbations will consist of all possible nonlinear gains 
lying in a sector [0, k}. In other words, we consider the configurations in Figures 51 and 52. 



Figure 51: 



Figure 52: 

The gain block <f> consists of all nonlinear time- varying gains <fr{t,a) satisfying 

</>(£, 0) = 0 for all t > 0 and 0 < a<t>(t,a) < ka 2 . 

This implies that <j>{t , a) is bounded by the lines <p = 0 and <j> = ka. Such nonlinearities are said to 
belong to a sector [0, k\. G(s) will be assumed to lie in the multilinear interval family G(s). The 
problem is to determine the largest size of the sector k for which robust stability is guaranteed. This 
is the multilinear version of the robust Lur’e problem. As before, the solution depends on the strict 
positive realness (SPR) properties of the family in question. First we have the following. 

Lemma 3. (Extremal SPR Properties) 

1) Let G(s) be the multilinear interval family defined in (58) and assume that G (s) is stable. Then 

inf inf Re(W(ju?)G(ju>)) = inf inf R e{W(ju)G(ju>)). 

GeG weR v ' v n Ge Ge^r 

2) IfC(s) is a controller that stabilizes the entire family G(s), then 

inf inf Re (Wtfw)C(jw^ = 

GeG 

inf inf Re (W(joj)C(ju)G(ju){l + C{ju))G(ju))~ l ) . 

G€Ge w€E 
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Theorem 21. (Nonlinear Stability margin) 

1) Let k* > 0 be defined by: 


sup/fc : - + inf inf Re(W (ju)G(ju))) > 0 1 
^ k G€Ge J 

then the closed loop system in Figure 51 is absolutely stable for all nonlinear gains <f> lying in 
the sector [0, A;*]. 

2) Let k* >0 be defined by: 

sup (fc : \ + inf inf Re (W(ju)C(juj)G(juj){l + C(jw)G(ju>)) _1 ) > ol 
J k GEGe wER J 

for the controller C, then the closed loop system in Figure 52 is absolutely stable for all non- 
linear gains (j) lying in the sector [0, k*]. 

Example 16. (Stability Sector for Multilinear Interval Systems) Consider the system used 
in Example 12. We plot the frequency domain image of the closed loop system transfer function 

^ 4S :i!|,|i(:,( 4 

Figure 53 shows that the nonlinear sector is given by £ = 20.5012. 



Figure 53: Frequency domain image of M(ju;) (Example 16) 


5.6 INTERVAL PLANTS AND DIAGONAL REPRESENTATION OF PERTURBA- 
TIONS 

In the robust control literature it has become customary to represent system perturbations in a 
signal flow diagram where the perturbations are “pulled” out and displayed in a feedback matrix 
A with independent diagonal or block- diagonal entries A* as shown below in Figure 54. M(s) is 
a suitably defined interconnection transfer function matrix of appropriate size and A is a block 
diagonal matrix containing all the perturbations affecting the system dynamics. In general this 
would include real parametric uncertainty as well as norm-bounded uncertainty blocks with the latter 
representing either actual unmodelled dynamics or fictitious performance blocks. The popularity of 
this representation is due to the fact that almost all types of uncertainties including parametric and 


51 



Figure 54: System with uncertainties represented in feedback form 


unstructured uncertainties can be accurately and explicitly represented in this framework. Moreover, 
the feedback representation allows us to use the Small Gain formulation and thus convert robust 
performance problems into robust stability problems. 

Although this formulation does not add anything to the the rich structure of interval systems 
already developed, it is nevertheless instructive to interpret the solution given by the Generalized 
Kharitonov Theorem in the linear and multilinear cases in this framework. We begin in the next 
subsection with the diagonal representation of a single interval system. This is subsequently extended 
to the case of multiple interval systems and multiple norm-bounded perturbation blocks. This type 
of mixed perturbation problem arises in systems with several performance specifications. We show 
how the boundary properties of the extremal systems Ge(s) established in the GKT play a role in 
the solution of these problems. 

5.6.1 Diagonal Feedback Representation of an Interval System 
Let G(s) be a interval plant 

G(s) = : (N(s) x D(s) € (N(s) x D(s))} (71) 

where N(s) and D(s) are interval polynomial families. The four Kharitonov polynomials associated 
with the interval polynomial D(s) can be written as follows: 

Kb(s) = D^n n (s) + J3SSW, KUs) = D™ n (s) + D**(s), 

K 3 d (s) = DZM + D^ n (s), K* d (s) = D*™(s) + D° m d « x (s). 

We know that, as far as frequency domain properties at s = ju> are concerned, the interval polynomial 
D(«) can be replaced by the following reduced 2-parameter family: 

D fl (s) = {£>(«) : D°(s) + XiD e (s) + X 2 D 0 (s) A* € [-1, 1]} (72) 

where 

DM = \ tfCSW - DZnis )) , DM = - - D^(s)) 

and 

D°(s) = \ [D™M + DZn n (s) + D°™M + D°^(s)] 

is the nominal polynomial. Similarly N(s) can be replaced by: 

N R (s) = {N(s) : N°(s) + A 3 N e (s) + A *N 0 (s) } € [-1, 1]} . (73) 

It is clear from the above that an interval plant is completely characterized in the frequency domain 
by the 4-parameter family Gr(s): 

G r (s ) = : (N(s) x D(s)) € (N r (s) x D r (s))| . (74) 

Now write y(s) = G(s)u(s), 


_ 7V«(s) + A 3 iVe(s) + A 4 iV 0 (s) 

V[) D°(s) + XMs) + x 2 dm 1 
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This may be rewritten as follows: 

y(s) = [D°(s )] _1 { ( N°(s ) + A 3 N e (s) + A Ms)) u(s ) - (A 1 D e (s) + A 2 D 0 (s)) y(s)} (75) 

or, suppressing s, as 

y = [(Z? 0 )" 1 ^ 0 + X 3 (D°)~ 1 N e + A 4 (r>°)- 1 ^ 0 ] u - [A^D 0 )" 1 ^ + A 2 (D°)- l D 0 ] y. (76) 

Corresponding to (76), we have the block diagram shown in Figure 55. 



Figure 55: Four parameter structure 


This input-output relation may be written as follows: 

■ -D e -D 0 N e N 0 N° 1 I" Wl ' 

-D e -D 0 N e N 0 N° w 2 

= (D 0 )- 1 0 0 0 0 D° w 3 (77) 

0 0 0 0 D° w 4 

_ -D e -D 0 N e No N° J [ 

with 

Wi = AjZi, i = 1,2,3, 4. (78) 

Let z := [z 1 ,z 2 ,Z3,z 4 ) t and w := [w\,w 2 ,w 3 , w 4 ] T . The above equations can then be rewritten as 
z(s) = M n (s)w(s) + M 12 (s)u(s), y(s) = M 2 i(s)«;(a) + M 22 (s)u(s), w(s) = Az(s). 

where 



M n (s): M 12 (s) 
M 2 i(s) : M 22 (s) 


M(s) = 

-\D°]~ l D e -[D°] -1 £> 0 [£>°] _1 JV e [D°] _1 iV 0 : [£>°] _1 ./V 0 
-\D°]- l D e -[D°] -1 D 0 [Z) 0 ] _1 i\T e [D°] _1 iV 0 : [D 0 ] _1 JV° 
0 0 0 0:1 
0 0 0 0:1 

~[D°}- l D e -[D^-'Do [L>°] _1 JV e [D 0 ] -1 AT o : [D 0 ]~ l N 0 


(79) 


(80) 


and 



Consequently, we have the configuration shown in Figure 56. 

The result established from GKT tells us that robust stability of the four parameter feedback 
structure given in Figure 56 is further equivalent to that of a reduced set of one-parameter extremal 
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Figure 56: 4-parameter feedback representation of an interval system 
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Figure 57: A typical element of the extremal set Ge(s) 


feedback structures in which each structure contains a single perturbation block. In Figure 57 we 
display a typical such extremal system structure where Ai is the only perturbation block and the 
remaining A’s are set to the vertex values -fl or -1. 

This is denoted in the equations by writing A 2 = A£, A 3 = A 3 and A 4 = AJ. The system equations 
assume the form 



M^s) M[ 2 {s) 
M^(s) M' 2 (s) 



Af'OO 


(81) 


with 


W\ — \\Z\ 


where 


M'n(s) = MM = (1 + A^D 0 ]- 1 ^)- 1 ^ 0 ]- 1 ^ 

M'i(s) = MM = (1 + [Z?°] - 1 Z) 0 ) - 1 [^>° ] “ 1 + WDT'Ne + A^D 0 ]- 1 ^. 

The interval plant G{s) can always be replaced by the above set of one-parameter structures 
Ge (s) regardless of the rest of the feedback system. If the intervals in which the parameters vary 
are not fixed apriort one can also determine their maximum permissible excursion using the same 
framework. The only difference in this case is that the Kharitonov polynomials, and therefore the 
system M' defined from them, have parameters which depend on the dilation parameter e as shown 
below in Figure 58. 
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Figure 58: An element of Ge(c, s) 


The limiting value e* of e is determined as the smallest value of e for which the control system 
containing one of the one parameter structures representing Ge(c, s) acquires a pole on the imaginary 
axis. Thus e* can be found using several methods which include the Segment Lemma and the 
Bounded Phase Conditions. Of course, by restricting attention to the class of controllers which 
satisfy the vertex conditions given in GKT we get further simplification in that all the A’s can be 
frozen at vertices, and e* can be found from the vertex systems. 
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The above discussion may be summarized by stating that if G(s) is an interval plant contained 
within a control system structure, we can always replace it with the set of one-parameter extremal 
systems Ge(s) for the purposes of determining worst case performance and stability margins and for 
carrying out worst case frequency response analysis. This observation holds also for linear interval 
systems and poly topic systems. In the case of multilinear interval systems we can again replace 
the plant set G(s) by the corresponding set of extremal systems Ge(s) which now consist of a set 
reduced dimensional multilinear systems in terms of the interval parameters A*. This can in turn 
be replaced by a polytopic family using the Mapping Theorem. Finally, an extremal one-parameter 
family can be constructed for this polytopic system as shown in this section. 

5.6.2 Interval Plant with Hoo Norm-Bounded Uncertainty 
Now let us consider the mixed uncertainty system in Figure 59. 



Figure 59: 4-parameter representation of interval system with unstructured uncertainty 

It consists of an interval plant with an unstructured perturbation block A lying in a Hoo ball of 
specified radius. As before, the interval plant can be represented in terms of four feedback parameters 
around a suitable interconnection transfer matrix M(s). We can now represent all the perturbations, 
parametric as well as unstructured, in feedback form around a suitably defined system P(s ). Thus, 
we have the structure shown in Figure 60 with an appropriate P(s). This four-parameter system 
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Figure 60: Feedback representation of mixed uncertainties 

can be reduced to the extremal set of single parameter uncertainty problems. A typical element of 
this extremal set is shown in Figure 61 with an appropriate P'(s). By eliminating the feedback loop 
associated with A this can also be represented as the structure shown in Figure 62. 

5.6.3 Multiple Interval Systems and Unstructured Blocks 
Let us now consider the general configuration given in Figure 63. 

M(s) is an interconnection transfer function matrix of suitable dimensions and A c/ (s) is a diagonal 
matrix containing all the unstructured system perturbations. The unstructured perturbations can 
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Figure 61: A typical element of the extremal set of systems 



Figure 63: System with uncertainties represented in feedback form 


be considered by introducing the class of perturbations D u (r) defined as 





> 

D» = < 

A y = 

< 

A r . 

: Aj € Ax 

4 


(82) 


where Ai is a ball in the space of i/oo real rational functions, denoted by RHoq , with radius 1 

A 1 = {A : A € RHoo, ||A||oo<l}. 

Parametric perturbations or uncertainty can be modeled by letting interval transfer functions rep- 
resent each physically distinct subsystem. These individual systems can then be “pulled out” and 
represented in feedback form as described earlier. The result of doing this is the general representa- 
tion of the system shown in Figure 64. 

In this representation, the A * 7 block represents multiple norm-bounded uncertainties, and the A p 
block accounts for parameter perturbations modeled as the set of independent interval systems G l (s) 
arranged in diagonal form: 


V G (q) 


{ 

' G i 


\ 


: Gi € GS 

\ 

i 

O' 

l 



(83) 


From the discussion of the last subsection we can immediately replace each interval system by a 
four-parameter set. Further we can also replace each interval system G l (s) by the corresponding 
extremal set G^: 


Dg(«) 


{ 

' Gx 

y 

\ 


: Gi G Gj; > 

\ 

G q . 

J 


(84) 
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Figure 64: Feedback representation of multiple interval plants and unstructured uncertainties 

For carrying out frequency domain analysis, the system in Figure 64 can therefore be replaced by 
Figure 65 where A UP belongs to the mixed complex-real perturbation class D UP defined as 



Figure 65: Feedback representation of mixed complex-real block diagonal uncertainties 

This representation can be used to formulate a variety of robust stability and performance problems. 
We illustrate this by considering two specific examples. 

Example 17. Consider the control system shown in Figure 66. 
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Figure 66: Robust performance problem (Example 17) 


In this system there are two independent subsystems containing parametric uncertainty which 
are modelled by the interval systems Gi{s) E G 1 ^) and (? 2 (s) € G 2 (s). Unstructured perturbations 
are represented by the block Ai , and appropriate scalings have been introduced to allow Ai to belong 
to Ai, the Hoo ball of radius unity. The specifications on the system performance are stated as the 
following requirements: 

1) Robust Stability The system must be stable for all Gi(s) E G x (s), G 2 (s) E G 2 (s), Ai E Ai 
and 


2) Robust Performance 
by unity: 


The worst case Hoc norm of the closed loop transfer function be bounded 
sup \\T ry (Gi , G 2 , Ai)||oo < 1, 

G l ,G 2 ,Ai 


where T ry is the transfer function from r to y. 

This problem can be cast in the framework developed above, and treated as a robust stability 
problem by introducing a fictitious performance block A 2 E Ai connected between r and y as 
shown in Figure 66. The equations describing the system in Figure 66 are 
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The “feedback” perturbations belong to the class 


D g ( 2) = { Gl 


g 2 


: G\ E G 1 ,G 2 E G‘ 


}■ 


and the interconnection matrix is 



M = 


0001-1 

1 0 0 0 0 

1 0 0 0 0 

0 110 0 

0 110 0 


Ai, A 2 6 Ai 


} 


Example 18. We consider the system in Figure 67. 

Parameter uncertainty is represented by the interval system Gi(s) G G 1 (s) and unstructured un- 
certainty consists of two independent blocks Ai,A 2 € Aj. Wi(s) and W 2 (s) represent frequency 
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Figure 67: Robust stability problem (Example 11.7) 


weightings and consist of fixed stable, proper minimum phase transfer functions. The system equa- 
tions are 

0 0 -1 -1 " w\ w\ G\ 0 0 0 zi 

0 0 —W\ —W\ w 2 w 2 _ 0 Ai 0 0 Z 2 

W 2 W 2 0 0 w 3 ’ w 3 “ 0 0 A 2 0 z 3 

1 1 0 0 J [ u J [ . u J ° 0 0 C J [ V m 

This configuration fits the general structure developed above with the “feedback” perturbation 
classes being 

D g (1) = {Gi : Gi 6 G 1 } and D G (2) = j [ Al ^ : A lf A 2 € j 

and the interconnection matrix 

“0 0 -1 -1 

0 0 -W, -Wi 

M ” W 2 W 2 0 0 

110 0 



5.6.4 Extremal Properties 

We have shown that parameter uncertainty can be modelled by interval systems and that both 
parametric and norm-bounded uncertainty can be represented in feedback form. Moreover, each 
interval system can be replaced, as far as worst case frequency domain analysis is concerned, by the 
corresponding reduced 4-parameter family. 

Let us now consider the system in Figure 64. Suppose that a stabilizing controller for the 
nominal system has been connected between y and u. The controller is stabilizing when the feedback 
perturbation is zero. Let M(s) denote the interconnection matrix with the controller attached. We 
now wish to determine if the control system is robustly stable under the feedback perturbations given. 
This is a real-complex mixed perturbation robust stability problem with the diagonal “feedback” 
perturbation matrix (see Figure 65) 



We make the standing assumption that the McMillan degree of the system remains invariant under 
the perturbations. The first observation we can make is obvious. 
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Theorem 22. The system in Figure 64 is Hurwitz stable for all D p 6 D G (m) and all D u G D G (r) 
if and only if the system in Figure 65 is Hurwitz stable for all D up G D UP . 

The above theorem actually applies more generally. For instance a more general formulation could 
allow for repeated interval system blocks. 

As a further simplification in the above result, we can replace each interval system G'(s) by 
the corresponding set of extremal systems G^(s). This corresponds to replacing each set of 4 A's in 
Figure 65 by their exposed edges. 

Theorem 23. The system in Figure 64 is Hurwitz stable for all D p G D G (m) and all D u G D c/ (r) 
if and only if the system in Figure 65 is Hurwitz stable for all Dp G Dg(m) and all D u G D G (r). 

Remark 8. The result stated above can be interpreted in terms of the worst case norm-bounded 
stability margin of the system over the interval parameter uncertainty set. For a prescribed D p , this 
stability margin is usually measured as a norm of the smallest destabilizing diagonal perturbation 
matrix A. The norm customarily used is the maximum singular value. The worst case value of this 
stability margin as D p ranges over the interval systems can be found by restricting D p to range 
over the extremal systems only. This is obviously a tremendous saving in computational effort. 
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6 INTERVAL MODEL IDENTIFICATION 

As demands for the control of flexible structures become more stringent, there arises an increasing 
need for more precise means of active vibration suppression. However, there is a growing need to 
develop controllers which are robust; that is the performance of the overall system needs to be 
insensitive to changes in the plant parameters. 

Although passive damping mechanisms such as in-house vibration absorbent materials are 
helpful, more stringent requirements call for active control schemes. Several types of actuator de- 
vices have been developed for applying forces to structures to meet this need. Results in the past 
have also shown that various design strategies exist for controlling vibration suppression actuators. 
Some successful local control designs using RMAs are described in various literatures. However, 
an important issue still remains open; that is maintaining a given level of performance in spite of 
changes in the plant dynamics. For the past 10 years, an extensive amount of research has been car- 
ried out in the area of parametric robust control, which deals with control systems with parametric 
uncertainties. 

Our aim is to apply these theoretical developments to a laboratory experiment, and to achieve 
meaningful robust control design and analysis of the vibration control of a flexible structure. There- 
fore, it can be guaranteed that the level of damping will be bounded inside a predicted range. The 
test structure is a scaled version of a typical flexible truss-like space structure and the parametric un- 
certainty is represented by mass added to various locations on the structure. This type of structural 
changes is common in many space-bound dynamic systems. For example, appendage articulation 
and decreasing fuel weight over the life of a satellite are possible system variations. While inertial 
changes are chosen for this investigation, this by no means limits the methodology presented herein 
to only inertial changes. Indeed, the effects of changes of any plant parameter such as structural 
damping or moduli can be accounted for the interval system techniques. 

The structural uncertainties can be modeled in different ways. In the approach, the 
system is described as a nominal model with its corresponding error bound transfer function. 
Minor structural changes may result in the changes of natural frequencies, damping ratios, and mode 
shapes. For instance, a small change of natural frequency of the lightly damping structures may 
result in significant frequency domain model error. Consequently, it may be necessary to find an 
Hqq error transfer function with large magnitude to bound the model error. Clearly, the choice of 
largely overbounded H 0 0 error bound transfer function results in overly conservative H ^ control 
design. 

A new model identification scheme has been introduced to overcome this difficulty. The tech- 
nique introduces a set of artificial parameters to fit the form of interval system framework so that a 
family of frequency response functions generated by this model contains the dynamics of the actual 
model. While this technique deals with a single frequency response function, in the work herein, we 
address the problem of obtaining an interval system model from a set of frequency response func- 
tions. These various frequency response functions can be attributed to parameter changes within 
the structure. Using this identified interval model, a robustness verification is performed to predict 
the open and closed loop uncertain system behavior. The interval system used in this study is given 
as a transfer function whose numerator and denominator are linear combinations of polynomials 
where the coefficients of the linear combinations are bounded with known values. 

The paper is organized as follows: first a description of our experimental system including 
the parameter uncertainty and control strategy is given. Next, a construction of the interval model 
which captures the system uncertainties due to change of mass is presented. Then detailed analysis 
of this interval model with the given controller is given. The results are compared with actual 
experimental data to validate the analysis. 

6.1 System Dynamics and Vibration Suppression 

6.1.1 Structural Dynamics 

A 10-bay aluminum truss structure shown in Figure. 68 is the test article. Each bay of the structure is 
a 0.5 meter cube made of hollow aluminum tubing joined at spherical nodes with threaded fasteners. 
The last three bays of the structure form a T-section. This is done so that the torsional modes of 
vibration will be more pronounced. The final structure has a mass of 31.5Ar^ and is cantilevered 
horizontally from a massive monolithic base. The ANSYS modeling package by Swanson Analysis is 
used to create a FEM model of the structure as shown in Figure. 69. A modal analysis is performed 
on the FEM model giving information about the structure’s modal frequencies and mode shapes. 
The results show that 7 modes are within the 100Hz bandwidth. 

To verify the FEM model, experimental frequency domain data are collected from the structure 
using a Tektronix 2642A Fourier Analyzer and compared with the analytically derived transfer 
function. The input to the structure consists of an impulse applied to Node 40 (see Figure. 69) of 
the structure, and the output is from an accelerometer mounted at Node 40 as well. The natural 
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Figure 68: Ten bay MeroForm truss with T-section. 


frequencies and damping ratios are obtained by using the Eigensystem Realization Algorithm (ERA) 
method. 

Table 1 shows the results for the first seven modes from the FEM as well as the identified 
models. From the FEM model we also extract a Guyan reduced order mass matrix M and stiffness 
matrix K for the structure. Hence the mathematical model of the structure including the ERA 
identified damping matrix D is represented by 

Mx + Dx + Kx = Bu. (86) 


6.1.2 Actuator Dynamics 

A reaction mass actuator is used to input force into the structure. This actuator uses a magnet and 
electrical-coil assembly to accelerate a mass. An amplified signal applied to the stationary electrical- 
coil creates a magnetic field. The constant field of the magnet inside the moving mass reacts against 
the new magnetic field causing force and hence movement of the mass. The mass rides on shafts 
which are limited to linear motion. The mass has a total usable travel length of about 2.5 cm. Any 
acceleration of the mass caused by force imparted from the coil results in an equal but opposite 
reaction force on the actuator’s housing. When the actuator is mounted on the structure, the force 
is imparted to the structure. Hence, the actuator is space-realizable, because it requires no inertial 
reference frame such as the earth. 

The actuator contains an internal non-contacting displacement sensor that helps to keep the 
mass centered with respect to its housing. The block diagram in Figure. 70 shows the local control 
scheme used to keep the mass centered. The PD type controller contains a proportional feedback 
loop, which acts as a spring, and a derivative feedback loop, which adds an equivalent viscous 
damping to the motion of the mass. The resulting system acts like a simple single-degree of freedom 
mass-spring-damper system. Notice a first-order filter is included to reduce the effects of high 
frequency noise in the damping loop. Realizing that the force output of the actuator is actually the 
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Figure 69: Finite element model of 10 bay truss including node number. 


Table 1: Results of FEM model and ERA identified model 


Mode 

FEM model 

ERA model 


Freq. 

Freq. 

Damping 


(Hz) 

(Hz) 

(0 

1 Vertical Bending 

9.22 

9.45 

.0120 

2 Horizontal Bending 

9.31 

9.70 

.0140 

3 Torsional 

18.13 

19.69 

.0080 

4 Vertical Bending 

47.43 

49.01 

.0020 

5 Horizontal Bending 

55.12 

60.30 

.0010 

6 Torsional 

84.73 

92.84 

.0006 

7 Mixed Motion 

87.48 

98.70 

.0003 


inertial acceleration of the mass, the output force of the actuator has the following dynamics. 

F(s) _ mrs 3 + ms 2 

V (s) mrs 3 + ms 2 + (c + kr)s + k 

where m is the mass of the reaction mass, k is the equivalent stiffness, c is the equivalent viscous 
damping constant, and r is the time constant of the noise filter. The transfer function represents 
the force output F(s) with respect to the voltage input V(s). 



Figure 70: Block diagram of the local centering control system of RMA. 


As shown in Figure. 68, the actuator is mounted between Node 39 and Node 40 of the structure. 
This location on the end of the structure is chosen because the modal participation factors of the 
first few modes of the structure are large at this location. 

It should be noted that several design considerations were enforced when building the actuator 
and setting the local feedback gains. First, the actuator is designed so that its mass is small compared 
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to the structure. The reaction mass is approximately 1.67 kg and the parasitic mass is l.6kg, so 
that the actuator is just over 10% of the structure’s mass. Note also that the efficiency ratio of the 
actuator itself is 51%. The damping ratio of the actuator affects the damping in the structure for 
both passively (without active control) and the active control design. This phenomenon has been 
studied and results in a compromise between damping in the structure and stability of the actuator. 
These considerations were taken into account when the internal damping in the actuator was chosen. 

6.1.3 Parametric Uncertainty 

The uncertainty of the structure we consider here is due to added mass to Node 44 and Node 17. We 
want to examine the structural changes that occur while each added mass varies from 0kg to 2.5kg. 
In order to observe the frequency domain behavior of the system with respect to the changes of these 
added masses, we take the following six samples of various added masses as shown in Table 2. 


Table 2: Masses added to the structure 


Case 

1 

2 

3 

4 

5 

6 

Node 17 

0 

0.5 

1.0 

1.5 

2.0 

2.5Kg 

Node 44 

0 

0.5 

1.0 

1.5 

2.0 

2.5Kg 


Our particular interest is the model uncertainty for the first horizontal mode and the first 
torsional mode with the natural frequencies ranging from 4 to 24Hz. The output measurement from 
an accelerometer located at Node 40 in the horizontal direction is used and the input is generated 
by applying the random signal to the RMA actuator. Figure. 71 illustrates the changes of the 
experimental data when these inertial parameters change. Therefore, it is reasonable to assume that 
for any value of added mass ranging from 0 to 2.5kg, the frequency response of the corresponding 
structure remains inside the thinnest envelope covering the frequency response of these six sampled 
structures. 



Figure 71: Experimental frequency response function with different added masses: From Right to 
Left 0 Kg, 1 Kg, 2 Kg, 3 Kg, ±Kg y 5 Kg. 

By combining Eqs. (86) and (87), we arrive at a transfer function expression of the actuator- 
structure system 

S 2 X(s,p) , , 

G{s ' v) = -T5T m 

where V^s) is the voltage across the actuator coils and s 2 X(s 1 p) is the response of an accelerometer. 
The unknown parameter, p, represents the uncertainty in the structural dynamics. Although we 
are using inertia to generate a family of perturbed systems, this by no means limits the techniques 
presented herein to inertial perturbations only. 
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6.1-4 Vibration Suppression Control Law 

The controller is designed using classical root loci techniques. However, it is later cast in terms of 
interval control systems to observe the robustness of the closed loop system. Vibration damping is 
accomplished with a local velocity feedback (LVF) control scheme. The velocity signal fed into a 
single actuator from an accelerometer location x at Node 40 horizontal direction is 

u = -K v x. (89) 

The control design is based on the model in Eq. (86). Here we consider the first horizontal bending 
and the first torsional modes. Feeding back velocity multiplied by a feedback gain, K v , the closed 
loop pole locations of the system are predicted. After varying the gain, a root locus plot is created. 
The resulting plot shows a trend of increased damping in all of the modes as the actuator gain is 
increased. 

To implement the controller, a negative velocity signal is fed to the actuator to “cancel out” 
the motion of the structure, i.e. reduce the velocity of the structure. An accelerometer is mounted on 
the housing of the actuator and connected to the structure. To get a velocity signal, the acceleration 
signal from the sensor must be integrated. However, integration of electrical signals results in 
exponential growth of any DC noise or DC offsets anywhere in the system. Therefore, we use a 
“velocity estimator” . The estimator is formulated as 

^ * , (90) 

o(s) S z + W C S + U)‘ 

where w c is the break frequency of the filter, V (a) is an estimated velocity signal and a(s) is the 
accelerometer input. The value, w c , used here is 0.5 Hz, and since the estimator acts like an ideal 
integrator at frequencies above 6w c , the feedback signal should be proportional to the velocity of 
the structure. 

Initially, the controllers are implemented on the structure with no added mass. The control is 
successful in increasing the system damping, but the question of robustness still remains unanswered. 
How much will slight changes in the system affect the closed loop control of the structure? The 
interval system model answers this question. 


6.2 Model Identification 


6.2.1 Model Structure 

The objective of this section is to construct a family of models such that its frequency response 
contains the six samples of experimental data shown in Figure. 71. Each data set consists of 900 
data points over the frequency ranging from 4 to 24 Hz. We first select the structure of the family 
of transfer functions as follow: 


G(s) = 



rcoOP + Eilt 1 a » r t( g ) 

do{s) + £”li 0i(Ji(s) 


€ [a ( ,at], E [S { , 



(91) 


where no ( 5 ) and do(s) are numerator and denominator of the nominal transfer function, and ri(s) 
and qi(s) are base polynomials of degree m and m - 1, respectively. The parameters ai and fii 
are the interval parameters which will be determined by their limiting values af and Clearly, 
upon the determination of the family G(s), we expect that the frequency response of every structure 
corresponding to each value of added mass ranging from 0 to 2.5kg remains inside that of the family 
G(s). 

In the next subsection, we give a brief description of how to determine the desired family 
of transfer functions G(s), equivalently determine a set of base polynomials (ri(s), <7i(s)), and the 
limiting values of interval parameters a* and /5» . This type of model is called an interval model. 


6.2.2 Interval Model Identification 

The algorithm to generate an interval model starts with identifying a model for each set of data. 
Any standard identification technique can be used for this purpose and here we use a least-squares 
based system identification technique which is widely used in practice. We first determine the order 
of the transfer function: 


9i(s) = 


ni(s) 

di(s) 


"y+Ci 8 ” 1 1 + •..+«} 
S m +di n _ 1 s m ~ 1 + ... +4 ’ 


* = 6 . 
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Once a judicious choice of the order of the transfer function is made, determination of the coefficients 
of each transfer function may be easily carried out by solving the least-squares problem. Let D, (w) 
be the * th sample data set consisting of 900 data points and W(u>) is a positive weighting function. 
The coefficients of the transfer function gi(s) are obtained by minimizing the following index: 

900 

Ji =^W(w k ) [R{u k ) 2 + J(wfc) 2 ] (92) 

k=\ 


where 

R(oJk) ■= Real [Di(uk)di(juj) - n<(jw)] , I(u k ) ■= Imag [Di{u k )di(juj) - rii{ju)] . 

We repeat this procedure for all six data sets to obtain the transfer functions 9i(s), i = 1, — ,6. A 
typical choice for the nominal model g 0 (s) is the average of all six of these identified models such as: 

Ho = K n i • ■ • n m] T and ^) = [<*0 d? • ■ • d^,_i] T (93) 

where 



The parameter perturbation vectors of the z th identified model are defined as 

An^ = nj — Ho a nd Ad { — d { ~ d^. (95) 

In order to determine the set of base polynomials r*(s) and qi(s) in Eq. (91), we employ a Singular 
Value Decomposition. To illustrate this SVD technique, we first form a matrix for the uncertainty 
of the denominator as 

AD = [Ad x A d 2 A4]. (96) 

Then we follow the following procedure: 

1) Compute the weighted uncertainty matrix AD W = W^ l AD where Wj is a diagonal matrix 

with its / th diagonal element being the standard deviation of the / th row of AD. 

2) Using SVD of A D w , compute the basis U of the perturbation matrix AD W (i.e., A D w = 

USV T ). 

3) Compute the basis for AD by Ud = WdU and let 

u d = [«* -• tXrfJ. (97) 

4) Compute the base polynomials corresponding to the denominator. 

m 

= E WO*'" 1 - * = 1,2, ,6 (98) 

1=1 

where y^. (/) is the / th component of the column vector u^.. 

5) Determine the coordinate vector of the uncertainty AD. 

[a £i ••• AgJ =U; 1 AD. (99) 

6) Determine the upper and lower bounds of the parameters. 

/?i + = max |Ag.(i), i — 1, 2, • • • , 6} , /?f = min {Ag.(J), i = 1,2, — , 6^ 

for l = 1, 2, • • ■ ,m. 

This interval model represents an optimal linear box of the determined polynomials to cover 
the perturbation. Similarly, we can obtain r*j(s), a“ and af for the numerator. 
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6.3 Interval Model Analysis 

The above procedure ensures that every transfer function <fr(s), i = 1, 2, • • • , 6 which represents the 
system with different added masses, is contained in the interval transfer function G(s). Consequently, 
our design and analysis will be performed with this interval model. However, the performance 
analysis of this type of transfer function cannot be done by traditional techniques such as Bode 
plots and the root locus method because the transfer function G{s) represents a family of systems 
rather than a single system. Here, we give a brief description of how to construct Bode and root 
locus plots for such systems. Since the interval transfer function represents a family of systems, its 
Bode plots will form envelopes and the root locus plot will be cluster plots. 

In this particular problem, the Bode magnitude plot is obtained as follows: First constructing 
the set of single parameter transfer functions, let us define the parameter vector which consists of 
interval parameters 


p := [tti * • • a m+ i 0\ • * * 0 m ] ’=[pi P2 • • • P2m+l] 
while the limiting values of each parameter are well-defined for 

Pi < Pi < pi, i = 1, 2, • ‘ • , 2m + 1. 

Then the interval model that represents the entire family may be regarded as G(s,p). A single 
parameter transfer function is obtained by fixing all parameters at each of their limiting intervals 
(called vertices ) except one parameter at a time. For example, 

[pi P 2 Pz ' * * ftm+ll) 


is a single parameter system. This leads us to a total of (2m + l)2 2m single parameter transfer 
functions. Using this set of single parameter systems, one can construct the Bode magnitude and 
phase envelopes as well as the plot of root clusters of the entire family G(s). 

Specifically, let V* be the set of vertices excluding such as 


Pi 1 P 2 > ’ ’ ' » Pi— 1 ’ Pi+l > * * * >P2m+l 
Pi > P 2 > * * * 1 Pi— 1 j Pi+l » ’ ’ * * P2m+1 


Vi 


= < 


[pi >pj > • * * >Pt— 1 1 Pi+l » * * * >P2m+l] » 

[pi, pi, •••, pi-i, Pi+i , j 


It is easy to see that Vi consists of 2 2m components. Let 

Pi = [Pi P2 ” * Pi-1 Pt+l * ‘ * P2m+l]* 

Then the boundaries of the Bode magnitude envelope are obtained by collecting 


sup sup |G (ju), [Xpi + (1 - A)pt >Pi e Vi])| 

i=3v,2m+l AC[0,1] 


Oliu . _ I 1 \ I 

hS*.4mi |g (i "' [Apr + (1 _ X)Pi ’ p< e Vil)l 

for every c 0 . The boundaries of the Bode phase envelope can also be obtained similarly. 

The problem of plotting the boundaries of root clusters is nothing but (2m + l)2 2m root locus 
problem. These boundaries coincide with those of the actual regions in which the poles of the family 
of systems lie. 


6.4 Experimental Results 

Utilizing the techniques described earlier, this section gives the design and analysis of the vibration 
suppression control of the ten-bay truss structure with added mass uncertainty. A reaction mass 
actuator (RM A) located between nodes 39 and 40 is used to excite and control the structure’s motion. 
An accelerometer located at Node 40 is used to measure the acceleration in the x direction at this 
position. Here we consider the model for the modes within the 4 — 24 Hz frequency range. There 
are two modes, one bending mode and one torsional mode, in this bandwidth. The experimental 
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frequency domain data are collected using a Tektronix 2642 A Fourier Analyzer. The 900 point 
frequency domain experimental data for the system with actuator dynamics is expressed as 

9M=9a(M)9i(uk), k = 1,2,..., 900 (100) 

where g t is the i th set of experimental data and gf represents the system dynamics excluding actuator 
dynamics, where the transfer function of the actuator is given by 

n a (s) 2.3853 x 10~ 3 s 3 + 1.6745s 2 

S ~ d a (s) ~ 2.3853 x 10~ 3 s 3 + 1.6745s 2 + 25.781s + 1251.3' 1 

Following the modeling procedure described, we first apply the least squares technique to each set 
of data g? to obtain the identified model for each case of added masses given in Table 2. Table 3 
shows the structural eigenvalues of these six identified models. 


Table 3: Eigenvalues of the identified models 


Model No. 

First mode 

Second mode 

1 

-1.3864±j61.094 

-1.4733±jl23.00 

2 

-1.3065±j59.930 

-1.3677±jll7.76 

3 

-1.1854±j58.892 

-1.4871±jll3.96 

4 

-1.1081±j57.902 

-1.4167±jll0.79 

5 

-1.0424±j56.929 

-1.4873±jl08.03 

6 

-0.9581±j55.724 

-1.5050±jl05.10 


After we obtain the identified structural models excluding actuator dynamics, we compute the 
nominal model as the average of these six models and we have 


n 0 {s) -3.463 x 10“ 3 s 4 - 1.173 x 10" 3 s 3 - 24.68s 2 

“ Ms) ~ s 4 + 5.241s 3 + 16255s 2 + 40225s + 4.412 x 10 7 ' 


( 102 ) 


We now apply the algorithm developed in Section 6.2 to obtain the structural interval model. Table 4 
shows the results of the interval model. 


Table 4: Results of interval model 


Qi 

<?2 

<73 

<74 


1.4326 x HT 1 

-2.5118 x lO" 1 1 

-1.8767 x W 2 - 

1.0399 x 10- 5 

s 2 8.6678 x 10 2 

5.5011 x 10 2 

-8.0878 x 10 2 

1.1301 x 10 3 

s 1 4.1568 x 10 3 

1.8545 x 10 3 

6.8649 x 10 3 

8.2205 x 10 2 

s° 4.0831 x 10 6 

2.5353 x 10 6 

-2.4265 x 10° 

6.1023 x 10 fl 

0+ 3.1117x10° 

2.8036 x 10- 1 

2.5123 x 10“ 2 

2.5140 x 10" 3 

0- -2.3396 x 10° 

-2.1055 x MT 1 - 

-3.3913 x 10" 2 - 

2.5161 x 10~ 3 

n 

T2 

rz 

^4 


-8.7360 x 10~ 4 

-3.2529 x 10"^ 

9.2234 x 10 

~ 5 0 

0 

s 3 -1.3370 x 10~ 4 

3.8897 x 10“ 4 

1.0547 x 10 

0 

lO 

1 

0 

s 2 -3.3372 x 10° 

-1.0513 x 10° 

-3.5316 x 1C 

1° 0 

0 

s 1 0 

0 

0 

1 

0 

s° 0 

0 

0 

0 

1 

a + 2.4339 x 10° 

1.0485 x 10° 

6.3217 x 10 

1 

to 

O 

0 

q- -1.8816x10° 

-1.3450 x 10° 

-3.5834 x IQ" 2 0 

0 


The results show that the uncertainty of the denominator is dominated by the uncertainty 
in the direction of the first singular vector. This phenomenon shows that the perturbation of the 
denominator coefficients is almost linear due to the changes of the added masses. The perturbation 
distributed in the direction of the fourth singular vector is around 1000 times smaller than that of the 
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first singular vector. It also shows that the model uncertainty of the numerator part is dominated by 
the uncertainty in the directions of the first two singular vectors. Figure. 72 shows the root clusters 
for the two structural modes of interest. Experimentally identified eigenvalues are also included for 
our study. These root clusters are generated by evaluating a family of transfer functions (an interval 
plant), and the method used here is based on the Edge theorem. This means that the poles of every 
system in the family are found in the regions bounded by the dotted lines in Figure. 72. 



Figure 72: Eigenvalue clusters of interval model and eigenvalues of experimental structure for various 
added masses. 

From the interval model of the structure, one can write the interval model of the system, which 
includes the actuator dynamics as follows: 

"<*(") [ n °( s ) + _ n a (s)n 0 (s) + S£t l a<n B (*)r < (*) 

&(S) ~ d 0 (s){do(3)+ZZlfa<(*)l daMdoW + LZlAdaWuM ' 

To verify the identified interval system by using the closed-loop experiment, we first design a 
local velocity feedback controller based on the root loci techniques for the vibration suppression of 
the structure. The transfer function of the controller is 

K( \ ~ - -110Kt>* 

S d k (s) s 2 + ns + 7r 2 

where K v is the velocity feedback gain. This controller is designed for the structure without added 
mass. This controller with K v = 30 increases the damping ratio of the first mode from 2.2% to 6.6% 
and the damping ratio of the second mode from 1.2% to 3.9%. The closed-loop interval system with 
the controller K(s) can be computed as 

dk{s) jn o (s)n 0 (s) + «i(na(s)r<(»))] 

d k (s) [d Q (s)do(s) + E’iLl 0i ( d a (*)*(«))] + T(s) 


where 


T(s) = n k (s) 


m+1 

n a (s)n 0 (s) + ^ an(n a (s)ri(s)) . 

x=i 


This interval system representation may be separated into the fixed and perturbed terms 


T(s) = 


t*(s) + 

/3(«) + EZi PM*) + ET=V <*Ms) ‘ 


(104) 


Figure. 73 depicts the magnitude envelope of the closed-loop interval system which shows that 
the envelope includes the magnitude plots of the experimental transfer functions. We also plot the 
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Figure 73: Magnitude envelope ( — ) of interval model and experimental data (..) for the closed-loop 
system frequency response functions. 


root clusters of the closed-loop interval system in Figure. 74. The damping for each mode of the 
closed-loop system is increased and it is verified by comparing the boundaries of the root clusters in 
Figures. 72 and 74. Thus, from this root cluster plot, one can predict the worst case damping of the 
closed-loop system, while the added mass varies anywhere in between 0 to 2.5 Kg. In other words, 
as long as the varied system remains inside the interval system family, the poles of the closed-loop 
system will remain inside the cluster regions shown in Figure. 74. 



Figure 74: Eigenvalue clusters of closed-loop interval model and eigenvalues of closed-loop experi- 
ment for K v = 60. 


6.5 Discussion and Remarks 

A vibration suppression control law was designed and implemented for an experimental truss struc- 
ture with parametric uncertainty. It is shown that the system identification technique developed 
here captures changes in the structural system response due to the parameter variations. Further, 
the identified interval model was used for robustness verification. That is, the performance bounds of 
the open and closed loop systems can be mapped with this identified interval model. The robustness 
verification is independent of the control design methodology and serves over the range of system 
behavior. 
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It is worthwhile to note that the number of single parameter systems required to construct 
the Bode and root cluster plots is relatively large. In fact, this number depends upon the number 
of interval polynomials for the Bode envelopes and the number of parameters for the root cluster 
plot. For example, if the interval transfer function is purely a ratio of two interval polynomials, 
only 32 systems are required to construct the Bode plots regardless of the number of parameters. 
The transfer function structure chosen here (see Eq. (91)) requires the same number of systems for 
both plots. This choice is made because our problem is low order and consists of multiple data sets. 
Nevertheless, since the vibration suppression is achieved by controlling the first few modes of the 
structure in most cases, the techniques are useful even if the structure is large. As a final note, 
extension to the case of multivariable system is still an open problem. 
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7 ROBUSTNESS ANALYSIS OF STATE SPACE MODEL 


Most of the results given in this book deal with polynomials containing parameter uncertainty. 
These results can be directly used when the system model is described by a transfer function whose 
coefficients contain the uncertain parameters. When the system model is described in the state space 
framework, the parameters appear as entries of the state space matrices. The polynomial theory 
can then be applied by first calculating the characteristic polynomial of the matrix as a function 
of its parameters. In this section, the aim is to provide some computational procedures which can 
determine robust stability and compute stability margins for the case in which parameters appear 
linearly in the state space matrices. Under the technical assumption that the perturbation structure 
is of unity rank the characteristic polynomial coefficients depend on the parameters in multilinear 
form. This allows us to use the Mapping Theorem to develop an effective computational technique 
to determine robust stability. This is illustrated with numerical examples. 

Next, we describe a Lyapunov based technique to handle perturbations of state space matrices. 
A stability region in parameter space is determined by this method. While this method does not 
require us to compute the characteristic polynomial, the stability region obtained from this method is 
conservative. However, the method gives a direct way to handle perturbations of state space matrices. 
Furthermore, with these formulas, the stability margin can be increased through optimization over 
the controller parameter space. This procedure is referred to as robustification. The results are 
compared with the previously described method, which used the Mapping Theorem, via examples. 

In the last part of the section we describe some formulas for the matrix stability margin for 
the real and complex cases and for some special classes of matrices. 

7.1 STATE SPACE PERTURBATIONS 

Consider the state space description of a linear system: 

x = Ax + Bu, y = Cx (105) 

with the output feedback control 

u = Ky. (106) 

The stability of the closed loop system is determined by the stability of the matrix M := A + BKC. 
We suppose that the matrices A, B, K, and C are subject to parameter perturbations. Let 

P ~ [Pi>P2, - - - , Pi] (1° 7 ) 

denote the parameters subject to uncertainty and set 

p = p° + Ap (108) 

where p° is the nominal parameter and Ap denotes a perturbation. Write 

M( p) = M( p° + Ap) = M(; p°) + AM(p°, Ap). 

Assuming that the entries of AM(p°, Ap) are linear functions of Ap, we can write 

AM(p°, Ap) = Api-Ei + A/>2-E2 H (- A piEi. (109) 

We shall say that the perturbation structure is of unity rank when each matrix Ei has unity 
rank. The special attraction of unity rank perturbation structures is the fact that in this case the 
coefficients of the characteristic polynomial of M are multilinear functions of Ap as shown below. 
When the Apj vary in intervals, this multilinear structure allows us to use the Mapping Theorem 
to develop an effective computational procedure to determine robust stability and stability margins 
in the parameter space p. The stability margin will be measured as the smallest norm of the 
vector Ap required to make M(p°, Ap) just unstable. 

In the robust stability literature, state space perturbations are often treated by supposing that 

AM = DUE (110) 

where the matrix U is regarded as a perturbation. In this formulation, one can calculate the smallest 
induced norm of U for which M + AM just becomes unstable. We remark that the parametric 
stability margin, defined as the vector norm of the smallest destabilizing vector Ap, has a physical 
significance in terms of the allowable perturbations of the parameter p. Such a direct significance 
cannot be attached to the matrix norm. Nevertheless, it has become customary to consider matrix 
valued perturbations and we accordingly define the matrix stability radius as the norm of the smallest 
destabilizing matrix. We give some formulas for the calculation of the matrix stability radius in the 
real and complex cases. 
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7.2 ROBUST STABILITY OF INTERVAL MATRICES 

We first establish that the unity rank perturbation structure leads to multilinear dependence of the 
characteristic polynomial coefficients on the parameter p. 

7.3 Unity Rank Perturbation Structure 

Let us suppose that 

M( p) := M( p°) +ApiEi + Ap 2 #2 H b (HI) 

Mo 


Lemma 4. Under the assumption that rank (Ei) = 1 for each i , the coefficients of the characteristic 
polynomial of M( p) are multilinear functions of p. 

7.3.1 Interval Matrix Stability via the Mapping Theorem 
The objectives are to solve the following problems: 

Problem 1 Determine if each matrix M{ p) remains stable as the parameter p ranges over given 
perturbation bounds p t r < Pi < pf, i = 1, • • • , / . 

Problem 2 With a stable Af(p°), determine the maximum value of e so that the matrix M(p) 
remains stable under all parameter perturbations ranging over p° { - Wit < Pi < 4- Wi€ for 

predetermined weights > 0. 

These problems may be effectively solved by using the fact that the characteristic polynomial 
of the matrix is a multilinear function of the parameters. 

The problem 1 can be solved by the following algorithm: 

Step 1: Determine the eigenvalues of the matrix M( p) with p fixed at each vertex of II. With 

this generate the characteristic polynomials corresponding to the vertices of II. 

Step 2: Verify the stability of the line segments connecting the vertex characteristic polynomials. 
This may be done by checking the Bounded Phase Condition or the Segment Lemma. 

We remark that the procedure outlined above does not require the determination of the characteristic 
polynomial as a function of the parameter p. It is enough to know that the function is multilinear. 
To determine the maximum value of e which solves the second problem, we may simply repeat the 
previous steps for incremental values of e. In fact, an upper bound l can be found as that value of 
e for which one of the vertices becomes just unstable. A lower bound e can be determined as the 
value of c for which a segment joining the vertices becomes unstable as follows: 

Step 1: Set e = e/2 

Step 2: Check the maximal phase differences of the vertex polynomials over the parameter box 

corresponding to e. 

Step 3: If the maximal phase difference is less than 7r radians, then increase e to e + (e - e)/2 for 

example, and repeat Step 2. 

Step 4: If the maximal phase difference is n radians or greater, then decrease e to c - (e - c)/2 

and repeat Step 2. 

Step 5: This iteration stops when the incremental step or decremental step becomes small enough. 

This gives a lower bound e and an upper bound e. 

If e and e are not close enough, we can refine the iteration by partitioning the interval uncer- 
tainty set into smaller boxes. 

The following examples illustrate this algorithm. 
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7.3.2 Numerical Examples 

Example 19. Consider the interval matrix: 



P° = \p\,P%,pI] = [-3, -2, 1] 

and 

Pi € [pf ,pf] = [-3 - e, -3 + c], P 2 € [pf ,p£] = [-2 - e, -2 + e], p 3 € [pf >P3 ] = t 1 “ £ - 1 + «)• 

The problem is to find the maximum value e* so that the matrix A(p) remains stable for all e 6 [0, €*]. 
Although the solution to this simple problem can be worked out analytically we work through the 
steps in some detail to illustrate the calculations involved. 

The characteristic polynomial of the matrix is: 

<5(s,p) = s(s-pi) -P2P3- 

In general this functional form is not required since only the vertex characteristic polynomials are 
needed and they can be found from the eigenvalues of the corresponding vertex matrices. Let us 
now compute the upper bound for e. We have eight vertex polynomials parametrized by e: 

A v (s) := {<S»(s,p) : P € V} 


where 

V := {{pi ,P 2 > P$ )) (Pi )p 2 ’Pi) 1 ( Pi >p 2 iPz)> ( Pi 1 P 2 1 P 3 )> (Pi ’ P 2 5^3 )’ 

(pf , pf 1 P 3 ) 1 (Pi 1 pf 1 pf ) 1 (pf > pf . pf ) } 

We found that the vertex polynomial 63 (s) has a jus root at e = 1. Thus we set e = 1. Using 
the multilinear version of GKT (Theorem 11.1), A (s) is robustly Hurwitz stable if and only if the 
following sets are Hurwitz stable: 

Lx = {a [A(a - pf) + (1 - A)(a -p+)] - pfpf : A € [0, 1]} 

L 2 = {a [A(a - pi) + (1 - A) (a -pf)] -P 2 P 3 '■ A € [Oil]} 

L 3 = {« [A(s -pf) + (1 - A)(s -pf)] -pjpf ■ A € [0,1]} 

U = {s [A(s - pf ) + (1 - A )(a - pf)] - pfpf : A € [0, 1 ]} 

M 5 = {a(a -pf ) - [A 1 P 2 + (1 — Ax)pJ] [A 2 pf + (1 - A 2 )p 3 + ] : (Ax, A 2 ) G [0, 1] x [0, 1]} 

M 6 = {a(a -pf) - [Axpf + (1 - Ax)pJ] [A 2 pf + (1 - A 2 )pJ] : (Ax, A 2 ) G [0, 1] x [0, 1]} . 

The Lj, i = 1,2, 3, 4 are line segments of polynomials, so we rewrite them as follows: 

Lx = A(a - pf - ptpt) + (1 - A)(a - pf ~ pfpf) 

U — A(a — Pi — P 2 P 3 " ) + (1 — A)(s — p+ — p 2 pf ) 

L 3 = A(a - pf -pfPz) + (1 - A )(s - pf -P 2 P 3 ) 

U = A(s -pf -pfpf ) + (1 — A) (a — pf -pfpf ) 

Now we need to generate the set of line segments that constructs the convex hull of the image sets 
of M 5 and M 6 . This can be done by connecting every pair of vertex polynomials. The vertex set 
corresponding to M 5 is: 

M 5 V (s) := {M 5 : (Ai , A 2 ) <E {( 0 , 0), (0, 1), (1, 0), (1, 1)}} . 

If we connect every pair of these vertex polynomials, we have the line segments: 

L 5 = a(a — Pi) - [A pfpf + (1 - A )pfpf] = A(s 2 -pfa - pfpf) + (1 - A) (a 2 - pfa - pfpf ) 

L 6 = a(a - pf) - [A pfpf + (1 - A)pf pf] = A(a 2 - pf a - pfpf ) + (1 - A)(a 2 - pf a - pf pf ) 

L-j = a(a - Pi) - [Apf pf + (1 - A)p 2 pf ] = A(a 2 - pf a - pfpf ) + (1 - A)(a 2 - pf s - p 2 "pf ) 

L 8 = a(a - Pi") - [A pfpf + (1 - A)pf pf ] = A(a 2 - pf a - pfpf ) + (1 - A)(a 2 - pf a - pfpf ) 

Lq = a(a — Pi) — [A pfpf + (1 - A)pf pf] = A(a 2 - pf a - pfpf ) + (1 - A)(a 2 - pf a - pfpf ) 

L 10 = s(a - pf ) - [A pfpf + (1 - A )pf pf ] = A(a 2 - pf a - pfpf) + (1 - A)(a 2 - pf a - pfpf). 
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Similarly, for M 6 we have the line segments: 

I'll = s(s-p+) - [Ap+p+ + (1 - A )p+pj] = A(s 2 -pis-ptpt) + (1 - A)(s 2 -p+s -P+P 3 -) 
L 12 = s(s - pi) - [A p+p+ + (1 - A)p 2 -p+] - A(s 2 - pi s-p+pi) + (1 - A)(s 2 - pf s - p^pj) 

L 13 = s(s - p+) - [Ap 2 P 3 + (1 - A)p 2 p 3 - ] = A(s 2 -pis- p 2 pj) + (1 - A)(a 2 - pfs - p 2 P 3 ) 

■I 14 = s (s — pi) — [Ap 2 p 3 + (1 — A)p 2 p 3 ] = A(s 2 — pi s — P 2 P 3 ) + (1 — A)(s — pi s — p 2 p 3 ) 

L 15 = s(s - pi) - [A p+pj 4- (1 - A)p 2 -p+] = A(s 2 - pfs - p^Pa ) + (1 - A)(s 2 -pis-p^pt) 

L 16 = s(s-pi)~ [A pi pt + (1 - A )p^pj] = A(s 2 - p^s - P^pt) + (1 - A)(s 2 - pis - p 2 p 3 "). 

The total number of line segments joining vertex pairs in the parameter space is 28. However the 
actual number of segments we checked is 16. This saving is due to the multilinear version of GKT 
(Theorem 11.1), which reduces the set to be checked. By testing these segments we find that they 
are all stable for e < 1. Since e = 1 corresponds to the instability of the vertex polynomial <S 3 (s) we 
conclude that the exact value of the stability margin is e = 1. The stability check of the segments 
can be carried out using the Segment Lemma. 

We can also solve this problem by using the Bounded Phase Condition. The set of vertices is 

A v («) = {Vi,i = 1 , 2 , • • • , 8 } 


where 


Vi (s) = s 2 - pi s - pipi = s 2 - 2s + 2, 
v 3 (s) = s 2 -pis -pi P3 =s 2 -2s, 
V>>(s) = s 2 - pf s - pipi = s 2 - 4s + 2, 
P 7 (s) = s 2 -Pie-pip^ = s 2 - 4s, 


V 2 (s) = s 2 - pi s - P 2 P 3 — s 2 - 2s + 6, 
V 4 {s) -s 2 -pis- P 2 P 3 = s 2 - 2s, 
Ve(s) = s 2 -Pis-p^pi = s 2 - 4s + 6, 
Vs(s) = s 2 -P 1 S-P 2 P 3 = s 2 - 4s. 


From the vertex set we see that the difference polynomials Vi(s) - Vj(s) are either constant, first 
order, antiHurwitz or of the form cs and each of these forms satisfy the conditions of the Vertex 
Lemma. Thus the stability of the vertices implies that of the edges. Thus the first encounter with 
instability can only occur on a vertex. This implies that the smallest value already found of e = 1 
for which a vertex becomes unstable is the correct value of the margin. 
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Figure 75: Maximum Phase Differences 4 >a v (w) ( in degrees) of Vertex Polynomials (Example 19) 
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This conclusion can also be verified by checking the phases of the vertices. Since there are 
duplicated vertices, we simply plot phase differences of six distinct vertices of A (ju) with c = 1. 
Figure 75 shows that the maximum phase difference plot as a function of frequency. This plot shows 
that the maximal phase difference never reaches 180 degrees confirming once again that e = 1 is 
indeed the true margin. 


Remark 9. The phase of a vertex which touches the origin cannot be determined and only the 
phase difference over the vertex set is meaningful. The phase condition can therefore only be used 
to determine whether the line segment excluding endpoints, intersects the origin. Thus, the stability 
of all the vertices must be verified independently. 


Example 20. Let 
^ = (A + BKC)x = 


-1 

0 
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' 1 

0 ' 
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0 

-2 

0 

+ 

0 

1 


0 
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0 

— 1 + ^2 


1 0 1 
0 1 0 
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where 


k\ G [ k j , k^ ] — [ €, e\ 


k 2 € [fc 2 ,*£] = [— e,e]. 


We first find all the vertex polynomials. 


A y(s) := {6i{s,ki) : k* € V, » = 1,2, 3, 4} 


V {(&!, ^ 2 ) : (fci" , & 2 O) {k\ , )> (^1 > ^2~)> (&i j ^2 )} • 

We found that the minimum value of c such that a vertex polynomial just becomes unstable is 1.75. 
Thus, l = 1.75. Then we proceed by checking either the phase condition or the Segment Lemma. 
If the Segment Lemma is applied, one must verify the stability of six segments joining the four 
vertices in A v (s). If the phase condition is applied, one must evaluate the phases of the four vertex 
polynomials and find the maximum phase difference at each frequency to observe whether it reaches 
180°. Again, the calculation shows that the smallest value of e that results in a segment becoming 
unstable is 1.75. Thus e - 1.75. This shows that the value obtained e = 1.75 is the true margin. 


The algorithm can also be applied to the robust stability problem for nonHurwitz regions. The 
following example illustrates the discrete time case. 

Example 2L Consider the discrete time system: 


x(k + 1) 


-0.5 0 k 2 

1 0.50 -1 

ki k\ 0.3 


x(k) 


For the nominal values of fc® = ko — 0, the system is Schur stable. We want to determine the 
maximum value of e* so that for alt parameters lying in the range 

&i € (— e*,e*) *2 € (-€*,€*) 

the system remains Schur stable. Using the procedure, we find the upper bound e = 0.2745 which 
is the minimum value of e which results in a vertex polynomial just becoming unstable. Figure 76 
shows that the maximum phase difference over all vertices at each 0 € [0, 2n) with e = e is less than 
180°. Thus we conclude from the Mapping Theorem that the exact parametric stability margin of 
this system is 0.2745. 
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